!>\file micro_mg2_0.F90
!! This file contains Morrison-Gettelman MP version 2.0 - update of MG
!! microphysics with prognostic precipitation.

!>\ingroup mg2mg3
!>\defgroup mg2_0_mp Morrison-Gettelman MP version 2.0
!! This module includes the MG microphysics version 2.0 - update of MG 
!! microphysics with prognostic precipitation.
!!
!!\author Andrew Gettelman, Hugh Morrison, Sean Santos
!! e-mail: morrison@ucar.edu, andrew@ucar.edu
!!\n Contributions from: Peter Caldwell, Xiaohong Liu and Steve Ghan
!!  
!! - Anning Cheng adopted for FV3GFS 9/29/2017
!! - Anning Cheng added GMAO ice conversion and Liu et al. Liquid water conversion 
!! in 10/12/2017
!! - S. Moorthi - Oct/Nov 2017 - optimized the code
!! - S. Moorthi - Nov 2017 - made the sedimentation quasi-implicit
!! - Version 2 history:
!!  - Sep 2011: Development begun
!!  - Feb 2013: Added of prognostic precipitation
!!  - Aug 2015: Published and released version (\cite Gettelman_2015_1 \cite Gettelman_2015_2 )
module micro_mg2_0
!---------------------------------------------------------------------------------
!
! NOTE: Modified to allow other microphysics packages (e.g. CARMA) to do ice
! microphysics in cooperation with the MG liquid microphysics. This is
! controlled by the do_cldice variable.
!
! If do_cldice is false, then MG microphysics should not update CLDICE or
! NUMICE; it is assumed that the other microphysics scheme will have updated
! CLDICE and NUMICE. The other microphysics should handle the following
! processes that would have been done by MG:
!   - Detrainment (liquid and ice)
!   - Homogeneous ice nucleation
!   - Heterogeneous ice nucleation
!   - Bergeron process
!   - Melting of ice
!   - Freezing of cloud drops
!   - Autoconversion (ice -> snow)
!   - Growth/Sublimation of ice
!   - Sedimentation of ice
!
! This option has not been updated since the introduction of prognostic
! precipitation, and probably should be adjusted to cover snow as well.
!
!---------------------------------------------------------------------------------
! Based on micro_mg (restructuring of former cldwat2m_micro)
! Author: Andrew Gettelman, Hugh Morrison.
! Contributions from: Xiaohong Liu and Steve Ghan
! December 2005-May 2010
! Description in: Morrison and Gettelman, 2008. J. Climate (MG2008)
!                 Gettelman et al., 2010 J. Geophys. Res. - Atmospheres (G2010)
! for questions contact Hugh Morrison, Andrew Gettelman
! e-mail: morrison@ucar.edu, andrew@ucar.edu
!---------------------------------------------------------------------------------
! Code comments added by HM, 093011
! General code structure:
!
! Code is divided into two main subroutines:
!   subroutine micro_mg_init --> initializes microphysics routine, should be called
!                                  once at start of simulation
!   subroutine micro_mg_tend --> main microphysics routine to be called each time step
!                                this also calls several smaller subroutines to calculate
!                                microphysical processes and other utilities
!
! List of external functions:
!   qsat_water --> for calculating saturation vapor pressure with respect to liquid water
!   qsat_ice --> for calculating saturation vapor pressure with respect to ice
!   gamma   --> standard mathematical gamma function
! .........................................................................
! List of inputs through use statement in fortran90:
! Variable Name                      Description                Units
! .........................................................................
! gravit          acceleration due to gravity                    m s-2
! rair            dry air gas constant for air                  J kg-1 K-1
! tmelt           temperature of melting point for water          K
! cpair           specific heat at constant pressure for dry air J kg-1 K-1
! rh2o            gas constant for water vapor                  J kg-1 K-1
! latvap          latent heat of vaporization                   J kg-1
! latice          latent heat of fusion                         J kg-1
! qsat_water      external function for calculating liquid water
!                 saturation vapor pressure/humidity              -
! qsat_ice        external function for calculating ice
!                 saturation vapor pressure/humidity              pa
! rhmini          relative humidity threshold parameter for
!                 nucleating ice                                  -
! .........................................................................
! NOTE: List of all inputs/outputs passed through the call/subroutine statement
!       for micro_mg_tend is given below at the start of subroutine micro_mg_tend.
!---------------------------------------------------------------------------------

! Procedures required:
! 1) An implementation of the gamma function (if not intrinsic).
! 2) saturation vapor pressure and specific humidity over water
! 3) svp over ice
use machine,        only : r8    => kind_phys
use funcphys,       only : fpvsl, fpvsi

!use wv_sat_methods, only: &
!    qsat_water => wv_sat_qsat_water, &
!    qsat_ice => wv_sat_qsat_ice

! Parameters from the utilities module.
use micro_mg_utils, only : pi,   omsm,  qsmall, mincld, rhosn, rhoi, &
                           rhow, rhows, ac,     bc,     ai,    bi,   &
                           aj,   bj,    ar,     br,     as,    bs,   &
                           mi0,  rising_factorial

implicit none
private
save

public :: micro_mg_init, micro_mg_tend, qcvar

! Switches for specification rather than prediction of droplet and crystal number
! note: number will be adjusted as needed to keep mean size within bounds,
! even when specified droplet or ice number is used
!
! If constant cloud ice number is set (nicons = .true.),
! then all microphysical processes except mass transfer due to ice nucleation
! (mnuccd) are based on the fixed cloud ice number. Calculation of
! mnuccd follows from the prognosed ice crystal number ni.

logical :: nccons  ! nccons = .true. to specify constant cloud droplet number
logical :: nicons  ! nicons = .true. to specify constant cloud ice number

! specified ice and droplet number concentrations
! note: these are local in-cloud values, not grid-mean
real(r8) :: ncnst  ! droplet num concentration when nccons=.true. (m-3)
real(r8) :: ninst  ! ice num concentration when nicons=.true. (m-3)

!=========================================================
! Private module parameters
!=========================================================

!Range of cloudsat reflectivities (dBz) for analytic simulator
real(r8), parameter :: csmin = -30._r8
real(r8), parameter :: csmax = 26._r8
real(r8), parameter :: mindbz = -99._r8
real(r8), parameter :: minrefl = 1.26e-10_r8    ! minrefl = 10._r8**(mindbz/10._r8)

! autoconversion size threshold for cloud ice to snow (m)
real(r8)            :: dcs, ts_au, ts_au_min, qcvar

! minimum mass of new crystal due to freezing of cloud droplets done
! externally (kg)
real(r8), parameter :: mi0l_min = 4._r8/3._r8*pi*rhow*(4.e-6_r8)**3
real(r8), parameter :: zero=0.0_r8, one=1.0_r8,  two=2.0_r8,  three=3.0_r8, &
                       four=4.0_r8, five=5.0_r8, six=6._r8,   half=0.5_r8,  &
                       ten=10.0_r8, forty=40.0_r8, oneo6=one/six

!=========================================================
! Constants set in initialization
!=========================================================

! Set using arguments to micro_mg_init
real(r8) :: g           ! gravity
real(r8) :: r           ! dry air gas constant
real(r8) :: rv          ! water vapor gas constant
real(r8) :: cpp         ! specific heat of dry air
real(r8) :: tmelt       ! freezing point of water (K)

! latent heats of:
real(r8) :: xxlv        ! vaporization
real(r8) :: xlf         ! freezing
real(r8) :: xxls        ! sublimation

real(r8) :: rhmini      ! Minimum rh for ice cloud fraction > 0.

! flags
logical  :: microp_uniform, do_cldice, use_hetfrz_classnuc

real(r8) :: rhosu       ! typical 850mn air density

real(r8) :: icenuct     ! ice nucleation temperature: currently -5 degrees C

real(r8) :: snowmelt    ! what temp to melt all snow: currently 2 degrees C
real(r8) :: rainfrze    ! what temp to freeze all rain: currently -5 degrees C

! additional constants to help speed up code
real(r8) :: gamma_br_plus1, gamma_bs_plus1, gamma_bi_plus1, gamma_bj_plus1
real(r8) :: gamma_br_plus4, gamma_bs_plus4, gamma_bi_plus4, gamma_bj_plus4
real(r8) :: xxlv_squared,   xxls_squared
real(r8) :: omeps, epsqs

character(len=16)  :: micro_mg_precip_frac_method  ! type of precipitation fraction method
real(r8)           :: micro_mg_berg_eff_factor     ! berg efficiency factor

logical  :: allow_sed_supersat ! Allow supersaturated conditions after sedimentation loop
logical  :: do_sb_physics ! do SB 2001 autoconversion or accretion physics
logical  :: do_ice_gmao
logical  :: do_liq_liu

!===============================================================================
contains
!===============================================================================

!>\ingroup mg2_0_mp
!! This subroutine calculates
subroutine micro_mg_init(                                         &
     kind, gravit, rair, rh2o, cpair, eps,                        &
     tmelt_in, latvap, latice,                                    &
     rhmini_in, micro_mg_dcs,ts_auto, mg_qcvar,                   &
     microp_uniform_in, do_cldice_in, use_hetfrz_classnuc_in,     &
     micro_mg_precip_frac_method_in, micro_mg_berg_eff_factor_in, &
     allow_sed_supersat_in, do_sb_physics_in,                     &
     do_ice_gmao_in, do_liq_liu_in,                               &
     nccons_in, nicons_in, ncnst_in, ninst_in)

  use micro_mg_utils, only : micro_mg_utils_init
  use wv_saturation,  only : gestbl

  !-----------------------------------------------------------------------
  !
  ! Purpose:
  ! initialize constants for MG microphysics
  !
  ! Author: Andrew Gettelman Dec 2005
  !
  !-----------------------------------------------------------------------

  integer,  intent(in)  :: kind         !< Kind used for reals
  real(r8), intent(in)  :: gravit
  real(r8), intent(in)  :: rair
  real(r8), intent(in)  :: rh2o
  real(r8), intent(in)  :: cpair
  real(r8), intent(in)  :: eps
! real(r8), intent(in)  :: fv
  real(r8), intent(in)  :: tmelt_in     !< Freezing point of water (K)
  real(r8), intent(in)  :: latvap
  real(r8), intent(in)  :: latice
  real(r8), intent(in)  :: rhmini_in    !< Minimum rh for ice cloud fraction > 0.
  real(r8), intent(in)  :: micro_mg_dcs
  real(r8), intent(in)  :: ts_auto(2)
  real(r8), intent(in)  :: mg_qcvar

  logical,  intent(in)  :: microp_uniform_in !< .true. = configure uniform for sub-columns
                                             !! .false. = use w/o sub-columns (standard)
  logical,  intent(in)  :: do_cldice_in      !< .true. = do all processes (standard)
                                             !! .false. = skip all processes affecting cloud ice
  logical,  intent(in)  :: use_hetfrz_classnuc_in !< use heterogeneous freezing

  character(len=16),intent(in)  :: micro_mg_precip_frac_method_in  !< type of precipitation fraction method
  real(r8),         intent(in)  :: micro_mg_berg_eff_factor_in     !< berg efficiency factor
  logical,  intent(in)  ::  allow_sed_supersat_in !< allow supersaturated conditions after sedimentation loop
  logical,  intent(in)  ::  do_sb_physics_in !< do SB autoconversion and accretion physics
  logical,  intent(in)  ::  do_ice_gmao_in
  logical,  intent(in)  ::  do_liq_liu_in

  logical,  intent(in)  :: nccons_in, nicons_in
  real(r8), intent(in)  :: ncnst_in,  ninst_in
  logical ip
  real(r8):: tmn, tmx, trice



  !-----------------------------------------------------------------------

  dcs       = micro_mg_dcs * 1.0e-6
  ts_au_min = ts_auto(1)
  ts_au     = ts_auto(2)
  qcvar     = mg_qcvar

  ! Initialize subordinate utilities module.
  call micro_mg_utils_init(kind, rair, rh2o, cpair, tmelt_in, latvap, latice, &
                           dcs)


  ! declarations for MG code (transforms variable names)

  g      = gravit              ! gravity
  r      = rair                ! dry air gas constant: note units(phys_constants are in J/K/kmol)
  rv     = rh2o                ! water vapor gas constant
  cpp    = cpair               ! specific heat of dry air
  tmelt  = tmelt_in
  rhmini = rhmini_in
  micro_mg_precip_frac_method = micro_mg_precip_frac_method_in
  micro_mg_berg_eff_factor    = micro_mg_berg_eff_factor_in
  allow_sed_supersat          = allow_sed_supersat_in
  do_sb_physics               = do_sb_physics_in
  do_ice_gmao                 = do_ice_gmao_in
  do_liq_liu                  = do_liq_liu_in

  nccons = nccons_in
  nicons = nicons_in
  ncnst  = ncnst_in
  ninst  = ninst_in

  ! latent heats

  xxlv = latvap         ! latent heat vaporization
  xlf  = latice         ! latent heat freezing
  xxls = xxlv + xlf     ! latent heat of sublimation

  ! flags
  microp_uniform      = microp_uniform_in
  do_cldice           = do_cldice_in
  use_hetfrz_classnuc = use_hetfrz_classnuc_in

  ! typical air density at 850 mb

  rhosu = 85000._r8 / (rair * tmelt)

  ! Maximum temperature at which snow is allowed to exist
  snowmelt = tmelt + two
  ! Minimum temperature at which rain is allowed to exist
  rainfrze = tmelt - forty

  ! Ice nucleation temperature
  icenuct  = tmelt - five

  ! Define constants to help speed up code (this limits calls to gamma function)
  gamma_br_plus1 = gamma(br+one)
  gamma_br_plus4 = gamma(br+four)
  gamma_bs_plus1 = gamma(bs+one)
  gamma_bs_plus4 = gamma(bs+four)
  gamma_bi_plus1 = gamma(bi+one)
  gamma_bi_plus4 = gamma(bi+four)
  gamma_bj_plus1 = gamma(bj+one)
  gamma_bj_plus4 = gamma(bj+four)

  xxlv_squared = xxlv * xxlv
  xxls_squared = xxls * xxls
  epsqs = eps
  omeps = one - epsqs
  tmn   = 173.16_r8
  tmx   = 375.16_r8
  trice = 35.00_r8
  ip    = .true.
  call gestbl(tmn ,tmx ,trice ,ip ,epsqs , latvap ,latice ,rh2o ,     &
              cpair ,tmelt_in )



end subroutine micro_mg_init

!===============================================================================
!microphysics routine for each timestep goes here...

!\ingroup mg2_0_mp
!> This subroutine is the main microphysics routine to be called each time step
!!
!! this also calls several smaller subroutines to calculate
!! microphysical processes and other utilities
subroutine micro_mg_tend (                                       &
     mgncol,             nlev,               deltatin,           &
     t,                            q,                            &
     qcn,                          qin,                          &
     ncn,                          nin,                          &
     qrn,                          qsn,                          &
     nrn,                          nsn,                          &
     relvar,                       accre_enhan_i,                &
     p,                            pdel,                         &
     cldn,    liqcldf,        icecldf,       qsatfac,            &
     qcsinksum_rate1ord,                                         &
     naai,                         npccnin,                      &
     rndst,                        nacon,                        &
     tlat,                         qvlat,                        &
     qctend,                       qitend,                       &
     nctend,                       nitend,                       &
     qrtend,                       qstend,                       &
     nrtend,                       nstend,                       &
     effc,               effc_fn,            effi,               &
     sadice,                       sadsnow,                      &
     prect,                        preci,                        &
     nevapr,                       evapsnow,                     &
     am_evp_st,                                                  &
     prain,                        prodsnow,                     &
     cmeout,                       deffi,                        &
     pgamrad,                      lamcrad,                      &
     qsout,                        dsout,                        &
     lflx,               iflx,                                   &
     rflx,               sflx,               qrout,              &
     reff_rain,                    reff_snow,                    &
     qcsevap,            qisevap,            qvres,              &
     cmeitot,            vtrmc,              vtrmi,              &
     umr,                          ums,                          &
     qcsedten,                     qisedten,                     &
     qrsedten,                     qssedten,                     &
     pratot,                       prctot,                       &
     mnuccctot,          mnuccttot,          msacwitot,          &
     psacwstot,          bergstot,           bergtot,            &
     melttot,                      homotot,                      &
     qcrestot,           prcitot,            praitot,            &
     qirestot,           mnuccrtot,          pracstot,           &
     meltsdttot,         frzrdttot,          mnuccdtot,          &
     nrout,                        nsout,                        &
     refl,               arefl,              areflz,             &
     frefl,              csrfl,              acsrfl,             &
     fcsrfl,                       rercld,                       &
     ncai,                         ncal,                         &
     qrout2,                       qsout2,                       &
     nrout2,                       nsout2,                       &
     drout2,                       dsout2,                       &
     freqs,                        freqr,                        &
     nfice,                        qcrat,                        &
     prer_evap, xlat, xlon, lprnt, iccn,  nlball)

  ! Constituent properties.
  use micro_mg_utils, only: mg_liq_props,  &
                            mg_ice_props,  &
                            mg_rain_props, &
                            mg_snow_props

  ! Size calculation functions.
  use micro_mg_utils, only: size_dist_param_liq,   &
                            size_dist_param_basic, &
                            avg_diameter

  ! Microphysical processes.
  use micro_mg_utils, only: ice_deposition_sublimation,    &
                            sb2001v2_liq_autoconversion,   &
                            sb2001v2_accre_cld_water_rain, &       
                            kk2000_liq_autoconversion,     &
                            ice_autoconversion,            &
                            immersion_freezing,            &
                            contact_freezing,              &
                            snow_self_aggregation,         &
                            accrete_cloud_water_snow,      &
                            secondary_ice_production,      &
                            accrete_rain_snow,             &
                            heterogeneous_rain_freezing,   &
                            accrete_cloud_water_rain,      &
                            self_collection_rain,          &
                            accrete_cloud_ice_snow,        &
                            evaporate_sublimate_precip,    &
                            bergeron_process_snow,         &
                            liu_liq_autoconversion,        &
                            gmao_ice_autoconversion,       &
                            size_dist_param_ice

  !Authors: Hugh Morrison, Andrew Gettelman, NCAR, Peter Caldwell, LLNL
  ! e-mail: morrison@ucar.edu, andrew@ucar.edu

  ! input arguments
  integer,  intent(in) :: mgncol                  ! number of microphysics columns
  integer,  intent(in) :: nlev                    ! number of layers
  integer,  intent(in) :: nlball(mgncol)          ! sedimentation start level
  real(r8), intent(in) :: xlat,xlon               ! number of layers
  real(r8), intent(in) :: deltatin                ! time step (s)
  real(r8), intent(in) :: t(mgncol,nlev)          ! input temperature (K)
  real(r8), intent(in) :: q(mgncol,nlev)          ! input h20 vapor mixing ratio (kg/kg)

  ! note: all input cloud variables are grid-averaged
  real(r8), intent(in) :: qcn(mgncol,nlev)        ! cloud water mixing ratio (kg/kg)
  real(r8), intent(in) :: qin(mgncol,nlev)        ! cloud ice mixing ratio (kg/kg)
  real(r8), intent(in) :: ncn(mgncol,nlev)        ! cloud water number conc (1/kg)
  real(r8), intent(in) :: nin(mgncol,nlev)        ! cloud ice number conc (1/kg)

  real(r8), intent(in) :: qrn(mgncol,nlev)        ! rain mixing ratio (kg/kg)
  real(r8), intent(in) :: qsn(mgncol,nlev)        ! snow mixing ratio (kg/kg)
  real(r8), intent(in) :: nrn(mgncol,nlev)        ! rain number conc (1/kg)
  real(r8), intent(in) :: nsn(mgncol,nlev)        ! snow number conc (1/kg)

  real(r8)             :: relvar(mgncol,nlev)     ! cloud water relative variance (-)
  real(r8)             :: accre_enhan(mgncol,nlev)! optional accretion
! real(r8), intent(in) :: relvar_i                ! cloud water relative variance (-)
  real(r8), intent(in) :: accre_enhan_i           ! optional accretion
                                                  ! enhancement factor (-)

  real(r8), intent(in) :: p(mgncol,nlev)          ! air pressure (pa)
  real(r8), intent(in) :: pdel(mgncol,nlev)       ! pressure difference across level (pa)

  real(r8), intent(in) :: cldn(mgncol,nlev)       ! cloud fraction (no units)
  real(r8), intent(in) :: liqcldf(mgncol,nlev)    ! liquid cloud fraction (no units)
  real(r8), intent(in) :: icecldf(mgncol,nlev)    ! ice cloud fraction (no units)
  real(r8), intent(in) :: qsatfac(mgncol,nlev)    ! subgrid cloud water saturation scaling factor (no units)
  logical, intent(in)  :: lprnt
  integer, intent(in)  :: iccn


  ! used for scavenging
  ! Inputs for aerosol activation
  real(r8), intent(inout) :: naai(mgncol,nlev)    ! ice nucleation number (from microp_aero_ts) (1/kg)
  real(r8), intent(in)    :: npccnin(mgncol,nlev) ! ccn activated number tendency (from microp_aero_ts) (1/kg*s)
! real(r8), intent(in)    :: npccn(mgncol,nlev)   ! ccn activated number tendency (from microp_aero_ts) (1/kg*s)
  real(r8)                :: npccn(mgncol,nlev)   ! ccn activated number tendency (from microp_aero_ts) (1/kg*s)

  ! Note that for these variables, the dust bin is assumed to be the last index.
  ! (For example, in CAM, the last dimension is always size 4.)
  real(r8), intent(in) :: rndst(mgncol,nlev,10)   ! radius of each dust bin, for contact freezing (from microp_aero_ts) (m)
  real(r8), intent(in) :: nacon(mgncol,nlev,10)   ! number in each dust bin, for contact freezing  (from microp_aero_ts) (1/m^3)
  
  ! output arguments

  real(r8), intent(out) :: qcsinksum_rate1ord(mgncol,nlev) ! 1st order rate for
  ! direct cw to precip conversion
  real(r8), intent(out) :: tlat(mgncol,nlev)      ! latent heating rate       (W/kg)
  real(r8), intent(out) :: qvlat(mgncol,nlev)     ! microphysical tendency qv (1/s)
  real(r8), intent(out) :: qctend(mgncol,nlev)    ! microphysical tendency qc (1/s)
  real(r8), intent(out) :: qitend(mgncol,nlev)    ! microphysical tendency qi (1/s)
  real(r8), intent(out) :: nctend(mgncol,nlev)    ! microphysical tendency nc (1/(kg*s))
  real(r8), intent(out) :: nitend(mgncol,nlev)    ! microphysical tendency ni (1/(kg*s))

  real(r8), intent(out) :: qrtend(mgncol,nlev)    ! microphysical tendency qr (1/s)
  real(r8), intent(out) :: qstend(mgncol,nlev)    ! microphysical tendency qs (1/s)
  real(r8), intent(out) :: nrtend(mgncol,nlev)    ! microphysical tendency nr (1/(kg*s))
  real(r8), intent(out) :: nstend(mgncol,nlev)    ! microphysical tendency ns (1/(kg*s))
  real(r8), intent(out) :: effc(mgncol,nlev)      ! droplet effective radius (micron)
  real(r8), intent(out) :: effc_fn(mgncol,nlev)   ! droplet effective radius, assuming nc = 1.e8 kg-1
  real(r8), intent(out) :: effi(mgncol,nlev)      ! cloud ice effective radius (micron)
  real(r8), intent(out) :: sadice(mgncol,nlev)    ! cloud ice surface area density (cm2/cm3)
  real(r8), intent(out) :: sadsnow(mgncol,nlev)   ! cloud snow surface area density (cm2/cm3)
  real(r8), intent(out) :: prect(mgncol)          ! surface precip rate (m/s)
  real(r8), intent(out) :: preci(mgncol)          ! cloud ice/snow precip rate (m/s)
  real(r8), intent(out) :: nevapr(mgncol,nlev)    ! evaporation rate of rain + snow (1/s)
  real(r8), intent(out) :: evapsnow(mgncol,nlev)  ! sublimation rate of snow (1/s)
  real(r8), intent(out) :: am_evp_st(mgncol,nlev) ! stratiform evaporation area (frac)
  real(r8), intent(out) :: prain(mgncol,nlev)     ! production of rain + snow (1/s)
  real(r8), intent(out) :: prodsnow(mgncol,nlev)  ! production of snow (1/s)
  real(r8), intent(out) :: cmeout(mgncol,nlev)    ! evap/sub of cloud (1/s)
  real(r8), intent(out) :: deffi(mgncol,nlev)     ! ice effective diameter for optics (radiation) (micron)
  real(r8), intent(out) :: pgamrad(mgncol,nlev)   ! ice gamma parameter for optics (radiation) (no units)
  real(r8), intent(out) :: lamcrad(mgncol,nlev)   ! slope of droplet distribution for optics (radiation) (1/m)
  real(r8), intent(out) :: qsout(mgncol,nlev)     ! snow mixing ratio (kg/kg)
  real(r8), intent(out) :: dsout(mgncol,nlev)     ! snow diameter (m)
  real(r8), intent(out) :: lflx(mgncol,2:nlev+1)  ! grid-box average liquid condensate flux (kg m^-2 s^-1)
  real(r8), intent(out) :: iflx(mgncol,2:nlev+1)  ! grid-box average ice condensate flux (kg m^-2 s^-1)
  real(r8), intent(out) :: rflx(mgncol,2:nlev+1)  ! grid-box average rain flux (kg m^-2 s^-1)
  real(r8), intent(out) :: sflx(mgncol,2:nlev+1)  ! grid-box average snow flux (kg m^-2 s^-1)
  real(r8), intent(out) :: qrout(mgncol,nlev)     ! grid-box average rain mixing ratio (kg/kg)
  real(r8), intent(out) :: reff_rain(mgncol,nlev) ! rain effective radius (micron)
  real(r8), intent(out) :: reff_snow(mgncol,nlev) ! snow effective radius (micron)
  real(r8), intent(out) :: qcsevap(mgncol,nlev)   ! cloud water evaporation due to sedimentation (1/s)
  real(r8), intent(out) :: qisevap(mgncol,nlev)   ! cloud ice sublimation due to sedimentation (1/s)
  real(r8), intent(out) :: qvres(mgncol,nlev)     ! residual condensation term to ensure RH < 100% (1/s)
  real(r8), intent(out) :: cmeitot(mgncol,nlev)   ! grid-mean cloud ice sub/dep (1/s)
  real(r8), intent(out) :: vtrmc(mgncol,nlev)     ! mass-weighted cloud water fallspeed (m/s)
  real(r8), intent(out) :: vtrmi(mgncol,nlev)     ! mass-weighted cloud ice fallspeed (m/s)
  real(r8), intent(out) :: umr(mgncol,nlev)       ! mass weighted rain fallspeed (m/s)
  real(r8), intent(out) :: ums(mgncol,nlev)       ! mass weighted snow fallspeed (m/s)
  real(r8), intent(out) :: qcsedten(mgncol,nlev)  ! qc sedimentation tendency (1/s)
  real(r8), intent(out) :: qisedten(mgncol,nlev)  ! qi sedimentation tendency (1/s)
  real(r8), intent(out) :: qrsedten(mgncol,nlev)  ! qr sedimentation tendency (1/s)
  real(r8), intent(out) :: qssedten(mgncol,nlev)  ! qs sedimentation tendency (1/s)

  ! microphysical process rates for output (mixing ratio tendencies) (all have units of 1/s)
  real(r8), intent(out) :: pratot(mgncol,nlev)    ! accretion of cloud by rain
  real(r8), intent(out) :: prctot(mgncol,nlev)    ! autoconversion of cloud to rain
  real(r8), intent(out) :: mnuccctot(mgncol,nlev) ! mixing ratio tend due to immersion freezing
  real(r8), intent(out) :: mnuccttot(mgncol,nlev) ! mixing ratio tend due to contact freezing
  real(r8), intent(out) :: msacwitot(mgncol,nlev) ! mixing ratio tend due to H-M splintering
  real(r8), intent(out) :: psacwstot(mgncol,nlev) ! collection of cloud water by snow
  real(r8), intent(out) :: bergstot(mgncol,nlev)  ! bergeron process on snow
  real(r8), intent(out) :: bergtot(mgncol,nlev)   ! bergeron process on cloud ice
  real(r8), intent(out) :: melttot(mgncol,nlev)   ! melting of cloud ice
  real(r8), intent(out) :: homotot(mgncol,nlev)   ! homogeneous freezing cloud water
  real(r8), intent(out) :: qcrestot(mgncol,nlev)  ! residual cloud condensation due to removal of excess supersat
  real(r8), intent(out) :: prcitot(mgncol,nlev)   ! autoconversion of cloud ice to snow
  real(r8), intent(out) :: praitot(mgncol,nlev)   ! accretion of cloud ice by snow
  real(r8), intent(out) :: qirestot(mgncol,nlev)  ! residual ice deposition due to removal of excess supersat
  real(r8), intent(out) :: mnuccrtot(mgncol,nlev) ! mixing ratio tendency due to heterogeneous freezing of rain to snow (1/s)
  real(r8), intent(out) :: pracstot(mgncol,nlev)  ! mixing ratio tendency due to accretion of rain by snow (1/s)
  real(r8), intent(out) :: meltsdttot(mgncol,nlev)! latent heating rate due to melting of snow  (W/kg)
  real(r8), intent(out) :: frzrdttot(mgncol,nlev) ! latent heating rate due to homogeneous freezing of rain (W/kg)
  real(r8), intent(out) :: mnuccdtot(mgncol,nlev) ! mass tendency from ice nucleation
  real(r8), intent(out) :: nrout(mgncol,nlev)     ! rain number concentration (1/m3)
  real(r8), intent(out) :: nsout(mgncol,nlev)     ! snow number concentration (1/m3)
  real(r8), intent(out) :: refl(mgncol,nlev)      ! analytic radar reflectivity
  real(r8), intent(out) :: arefl(mgncol,nlev)     ! average reflectivity will zero points outside valid range
  real(r8), intent(out) :: areflz(mgncol,nlev)    ! average reflectivity in z.
  real(r8), intent(out) :: frefl(mgncol,nlev)     ! fractional occurrence of radar reflectivity
  real(r8), intent(out) :: csrfl(mgncol,nlev)     ! cloudsat reflectivity
  real(r8), intent(out) :: acsrfl(mgncol,nlev)    ! cloudsat average
  real(r8), intent(out) :: fcsrfl(mgncol,nlev)    ! cloudsat fractional occurrence of radar reflectivity
  real(r8), intent(out) :: rercld(mgncol,nlev)    ! effective radius calculation for rain + cloud
  real(r8), intent(out) :: ncai(mgncol,nlev)      ! output number conc of ice nuclei available (1/m3)
  real(r8), intent(out) :: ncal(mgncol,nlev)      ! output number conc of CCN (1/m3)
  real(r8), intent(out) :: qrout2(mgncol,nlev)    ! copy of qrout as used to compute drout2
  real(r8), intent(out) :: qsout2(mgncol,nlev)    ! copy of qsout as used to compute dsout2
  real(r8), intent(out) :: nrout2(mgncol,nlev)    ! copy of nrout as used to compute drout2
  real(r8), intent(out) :: nsout2(mgncol,nlev)    ! copy of nsout as used to compute dsout2
  real(r8), intent(out) :: drout2(mgncol,nlev)    ! mean rain particle diameter (m)
  real(r8), intent(out) :: dsout2(mgncol,nlev)    ! mean snow particle diameter (m)
  real(r8), intent(out) :: freqs(mgncol,nlev)     ! fractional occurrence of snow
  real(r8), intent(out) :: freqr(mgncol,nlev)     ! fractional occurrence of rain
  real(r8), intent(out) :: nfice(mgncol,nlev)     ! fractional occurrence of ice
  real(r8), intent(out) :: qcrat(mgncol,nlev)     ! limiter for qc process rates (1=no limit --> 0. no qc)

  real(r8), intent(out) :: prer_evap(mgncol,nlev)


  ! Tendencies calculated by external schemes that can replace MG's native
  ! process tendencies.

  ! Used with CARMA cirrus microphysics
  ! (or similar external microphysics model)
  ! real(r8), intent(in) :: tnd_qsnow(:,:)        ! snow mass tendency (kg/kg/s)
  ! real(r8), intent(in) :: tnd_nsnow(:,:)        ! snow number tendency (#/kg/s)
  ! real(r8), intent(in) :: re_ice(:,:)           ! ice effective radius (m)

  ! From external ice nucleation.
  !real(r8), intent(in) :: frzimm(:,:)            ! Number tendency due to immersion freezing (1/cm3)
  !real(r8), intent(in) :: frzcnt(:,:)            ! Number tendency due to contact freezing (1/cm3)
  !real(r8), intent(in) :: frzdep(:,:)            ! Number tendency due to deposition nucleation (1/cm3)

  ! local workspace
  ! all units mks unless otherwise stated

  ! local copies of input variables
  real(r8) :: qc(mgncol,nlev)         ! cloud liquid mixing ratio (kg/kg)
  real(r8) :: qi(mgncol,nlev)         ! cloud ice mixing ratio (kg/kg)
  real(r8) :: nc(mgncol,nlev)         ! cloud liquid number concentration (1/kg)
  real(r8) :: ni(mgncol,nlev)         ! cloud liquid number concentration (1/kg)
  real(r8) :: qr(mgncol,nlev)         ! rain mixing ratio (kg/kg)
  real(r8) :: qs(mgncol,nlev)         ! snow mixing ratio (kg/kg)
  real(r8) :: nr(mgncol,nlev)         ! rain number concentration (1/kg)
  real(r8) :: ns(mgncol,nlev)         ! snow number concentration (1/kg)

  ! general purpose variables
  real(r8) :: deltat                  ! sub-time step (s)
  real(r8) :: oneodt                  ! one / deltat
  real(r8) :: mtime                   ! the assumed ice nucleation timescale

  ! physical properties of the air at a given point
  real(r8) :: rho(mgncol,nlev)        ! density (kg m-3)
  real(r8) :: rhoinv(mgncol,nlev)     ! one / density (kg m-3)
  real(r8) :: dv(mgncol,nlev)         ! diffusivity of water vapor
  real(r8) :: mu(mgncol,nlev)         ! viscosity
  real(r8) :: sc(mgncol,nlev)         ! schmidt number
  real(r8) :: rhof(mgncol,nlev)       ! density correction factor for fallspeed

  ! cloud fractions
  real(r8) :: precip_frac(mgncol,nlev)! precip fraction assuming maximum overlap
  real(r8) :: cldm(mgncol,nlev)       ! cloud fraction
  real(r8) :: icldm(mgncol,nlev)      ! ice cloud fraction
  real(r8) :: lcldm(mgncol,nlev)      ! liq cloud fraction
  real(r8) :: qsfm(mgncol,nlev)       ! subgrid cloud water saturation scaling factor

  ! mass mixing ratios
  real(r8) :: qcic(mgncol,nlev)       ! in-cloud cloud liquid
  real(r8) :: qiic(mgncol,nlev)       ! in-cloud cloud ice
  real(r8) :: qsic(mgncol,nlev)       ! in-precip snow
  real(r8) :: qric(mgncol,nlev)       ! in-precip rain

  ! number concentrations
  real(r8) :: ncic(mgncol,nlev)       ! in-cloud droplet
  real(r8) :: niic(mgncol,nlev)       ! in-cloud cloud ice
  real(r8) :: nsic(mgncol,nlev)       ! in-precip snow
  real(r8) :: nric(mgncol,nlev)       ! in-precip rain
  ! maximum allowed ni value
  real(r8) :: nimax(mgncol,nlev)

  ! Size distribution parameters for:
  ! cloud ice
  real(r8) :: lami(mgncol,nlev)       ! slope
  real(r8) :: n0i(mgncol,nlev)        ! intercept
  ! cloud liquid
  real(r8) :: lamc(mgncol,nlev)       ! slope
  real(r8) :: pgam(mgncol,nlev)       ! spectral width parameter
  ! snow
  real(r8) :: lams(mgncol,nlev)       ! slope
  real(r8) :: n0s(mgncol,nlev)        ! intercept
  ! rain
  real(r8) :: lamr(mgncol,nlev)       ! slope
  real(r8) :: n0r(mgncol,nlev)        ! intercept

  ! Rates/tendencies due to:

  ! Instantaneous snow melting
  real(r8) :: minstsm(mgncol,nlev)    ! mass mixing ratio
  real(r8) :: ninstsm(mgncol,nlev)    ! number concentration
  ! Instantaneous rain freezing
  real(r8) :: minstrf(mgncol,nlev)    ! mass mixing ratio
  real(r8) :: ninstrf(mgncol,nlev)    ! number concentration

  ! deposition of cloud ice
  real(r8) :: vap_dep(mgncol,nlev)    ! deposition from vapor to ice PMC 12/3/12
  ! sublimation of cloud ice
  real(r8) :: ice_sublim(mgncol,nlev) ! sublimation from ice to vapor PMC 12/3/12
  ! ice nucleation
  real(r8) :: nnuccd(mgncol,nlev)     ! number rate from deposition/cond.-freezing
  real(r8) :: mnuccd(mgncol,nlev)     ! mass mixing ratio
  ! freezing of cloud water
  real(r8) :: mnuccc(mgncol,nlev)     ! mass mixing ratio
  real(r8) :: nnuccc(mgncol,nlev)     ! number concentration
  ! contact freezing of cloud water
  real(r8) :: mnucct(mgncol,nlev)     ! mass mixing ratio
  real(r8) :: nnucct(mgncol,nlev)     ! number concentration
  ! deposition nucleation in mixed-phase clouds (from external scheme)
  real(r8) :: mnudep(mgncol,nlev)     ! mass mixing ratio
  real(r8) :: nnudep(mgncol,nlev)     ! number concentration
  ! ice multiplication
  real(r8) :: msacwi(mgncol,nlev)     ! mass mixing ratio
  real(r8) :: nsacwi(mgncol,nlev)     ! number concentration
  ! autoconversion of cloud droplets
  real(r8) :: prc(mgncol,nlev)        ! mass mixing ratio
  real(r8) :: nprc(mgncol,nlev)       ! number concentration (rain)
  real(r8) :: nprc1(mgncol,nlev)      ! number concentration (cloud droplets)
  ! self-aggregation of snow
  real(r8) :: nsagg(mgncol,nlev)      ! number concentration
  ! self-collection of rain
  real(r8) :: nragg(mgncol,nlev)      ! number concentration
  ! collection of droplets by snow
  real(r8) :: psacws(mgncol,nlev)     ! mass mixing ratio
  real(r8) :: npsacws(mgncol,nlev)    ! number concentration
  ! collection of rain by snow
  real(r8) :: pracs(mgncol,nlev)      ! mass mixing ratio
  real(r8) :: npracs(mgncol,nlev)     ! number concentration
  ! freezing of rain
  real(r8) :: mnuccr(mgncol,nlev)     ! mass mixing ratio
  real(r8) :: nnuccr(mgncol,nlev)     ! number concentration
  ! freezing of rain to form ice (mg add 4/26/13)
  real(r8) :: mnuccri(mgncol,nlev)    ! mass mixing ratio
  real(r8) :: nnuccri(mgncol,nlev)    ! number concentration
  ! accretion of droplets by rain
  real(r8) :: pra(mgncol,nlev)        ! mass mixing ratio
  real(r8) :: npra(mgncol,nlev)       ! number concentration
  ! autoconversion of cloud ice to snow
  real(r8) :: prci(mgncol,nlev)       ! mass mixing ratio
  real(r8) :: nprci(mgncol,nlev)      ! number concentration
  ! accretion of cloud ice by snow
  real(r8) :: prai(mgncol,nlev)       ! mass mixing ratio
  real(r8) :: nprai(mgncol,nlev)      ! number concentration
  ! evaporation of rain
  real(r8) :: pre(mgncol,nlev)        ! mass mixing ratio
  ! sublimation of snow
  real(r8) :: prds(mgncol,nlev)       ! mass mixing ratio
  ! number evaporation
  real(r8) :: nsubi(mgncol,nlev)      ! cloud ice
  real(r8) :: nsubc(mgncol,nlev)      ! droplet
  real(r8) :: nsubs(mgncol,nlev)      ! snow
  real(r8) :: nsubr(mgncol,nlev)      ! rain
  ! bergeron process
  real(r8) :: berg(mgncol,nlev)       ! mass mixing ratio (cloud ice)
  real(r8) :: bergs(mgncol,nlev)      ! mass mixing ratio (snow)


  ! fallspeeds
  ! number-weighted
  real(r8) :: uns(mgncol,nlev)        ! snow
  real(r8) :: unr(mgncol,nlev)        ! rain
  ! air density corrected fallspeed parameters
  real(r8) :: arn(mgncol,nlev)        ! rain
  real(r8) :: asn(mgncol,nlev)        ! snow
  real(r8) :: acn(mgncol,nlev)        ! cloud droplet
  real(r8) :: ain(mgncol,nlev)        ! cloud ice
  real(r8) :: ajn(mgncol,nlev)        ! cloud small ice

  ! Mass of liquid droplets used with external heterogeneous freezing.
  real(r8) :: mi0l(mgncol)

  ! saturation vapor pressures
  real(r8) :: esl(mgncol,nlev)        ! liquid
  real(r8) :: esi(mgncol,nlev)        ! ice
  real(r8) :: esn                     ! checking for RH after rain evap

  ! saturation vapor mixing ratios
  real(r8) :: qvl(mgncol,nlev)        ! liquid
  real(r8) :: qvi(mgncol,nlev)        ! ice
  real(r8) :: qvn                     ! checking for RH after rain evap

  ! relative humidity
  real(r8) :: relhum(mgncol,nlev)

  ! parameters for cloud water and cloud ice sedimentation calculations
  real(r8) :: fc(mgncol,nlev)
  real(r8) :: fnc(mgncol,nlev)
  real(r8) :: fi(mgncol,nlev)
  real(r8) :: fni(mgncol,nlev)

  real(r8) :: fr(mgncol,nlev)
  real(r8) :: fnr(mgncol,nlev)
  real(r8) :: fs(mgncol,nlev)
  real(r8) :: fns(mgncol,nlev)

  real(r8) :: faloutc(nlev)
  real(r8) :: faloutnc(nlev)
  real(r8) :: falouti(nlev)
  real(r8) :: faloutni(nlev)

  real(r8) :: faloutr(nlev)
  real(r8) :: faloutnr(nlev)
  real(r8) :: falouts(nlev)
  real(r8) :: faloutns(nlev)

  real(r8) :: faltndc
  real(r8) :: faltndnc
  real(r8) :: faltndi
  real(r8) :: faltndni
  real(r8) :: faltndqie
  real(r8) :: faltndqce

  real(r8) :: faltndr
  real(r8) :: faltndnr
  real(r8) :: faltnds
  real(r8) :: faltndns

  real(r8) :: rainrt(mgncol,nlev)     ! rain rate for reflectivity calculation

  ! dummy variables
  real(r8) :: dum
  real(r8) :: dum1
  real(r8) :: dum2
  real(r8) :: dum3
  real(r8) :: dumni0
  real(r8) :: dumns0
  real(r8) :: tx1, tx2, tx3, tx4, tx5, tx6, tx7, grho
  ! dummies for checking RH
  real(r8) :: qtmp
  real(r8) :: ttmp
  ! dummies for conservation check
  real(r8) :: ratio
  real(r8) :: tmpfrz
  ! dummies for in-cloud variables
  real(r8) :: dumc(mgncol,nlev)   ! qc
  real(r8) :: dumnc(mgncol,nlev)  ! nc
  real(r8) :: dumi(mgncol,nlev)   ! qi
  real(r8) :: dumni(mgncol,nlev)  ! ni
  real(r8) :: dumr(mgncol,nlev)   ! rain mixing ratio
  real(r8) :: dumnr(mgncol,nlev)  ! rain number concentration
  real(r8) :: dums(mgncol,nlev)   ! snow mixing ratio
  real(r8) :: dumns(mgncol,nlev)  ! snow number concentration
  ! Array dummy variable
 !real(r8) :: dum_2D(mgncol,nlev)
  real(r8) :: pdel_inv(mgncol,nlev)
  real(r8) :: ts_au_loc(mgncol)

  ! loop array variables
  ! "i" and "k" are column/level iterators for internal (MG) variables
  ! "n" is used for other looping (currently just sedimentation)
  integer i, k, n

  ! number of sub-steps for loops over "n" (for sedimentation)
  integer nstep, mdust, nlb, nstep_def

  ! Varaibles to scale fall velocity between small and regular ice regimes.
  real(r8) :: irad, ifrac, tsfac
! logical, parameter  :: do_ice_gmao=.false., do_liq_liu=.false.
! logical, parameter  :: do_ice_gmao=.true., do_liq_liu=.true.
! logical, parameter  :: do_ice_gmao=.true., do_liq_liu=.false.
  real(r8), parameter :: qimax=0.010, qimin=0.005, qiinv=one/(qimax-qimin)
!                        ts_au_min=180.0

  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


  ! Process inputs

  ! assign variable deltat to deltatin
  deltat    = deltatin
  oneodt    = one / deltat
! nstep_def = max(1, nint(deltat/20))
  nstep_def = max(1, nint(deltat/5))
! tsfac     = log(ts_au/ts_au_min) * qiinv

  ! Copies of input concentrations that may be changed internally.
  do k=1,nlev
    do i=1,mgncol
      qc(i,k) = qcn(i,k)
      nc(i,k) = ncn(i,k)
      qi(i,k) = qin(i,k)
      ni(i,k) = nin(i,k)
      qr(i,k) = qrn(i,k)
      nr(i,k) = nrn(i,k)
      qs(i,k) = qsn(i,k)
      ns(i,k) = nsn(i,k)
    enddo
  enddo

  ! cldn: used to set cldm, unused for subcolumns
  ! liqcldf: used to set lcldm, unused for subcolumns
  ! icecldf: used to set icldm, unused for subcolumns

  if (microp_uniform) then
     ! subcolumns, set cloud fraction variables to one
     ! if cloud water or ice is present, if not present
     ! set to mincld (mincld used instead of zero, to prevent
     ! possible division by zero errors).

    do k=1,nlev
      do i=1,mgncol

        if (qc(i,k) >= qsmall) then
          lcldm(i,k) = one
        else
          lcldm(i,k) = mincld
        endif

        if (qi(i,k) >= qsmall) then
          icldm(i,k) = one
        else
          icldm(i,k) = mincld
        endif

        cldm(i,k) = max(icldm(i,k), lcldm(i,k))
!       qsfm(i,k) = one
        qsfm(i,k) = qsatfac(i,k)
      enddo
    enddo

  else         ! get cloud fraction, check for minimum
    do k=1,nlev
      do i=1,mgncol
        cldm(i,k)  = max(cldn(i,k), mincld)
        lcldm(i,k) = max(liqcldf(i,k), mincld)
        icldm(i,k) = max(icecldf(i,k), mincld)
        qsfm(i,k)  = qsatfac(i,k)
      enddo
    enddo
  end if

! if (lprnt) write(0,*)' cldm=',cldm(1,nlev-20:nlev)
! if (lprnt) write(0,*)' liqcldf=',liqcldf(1,nlev-20:nlev)
! if (lprnt) write(0,*)' lcldm=',lcldm(1,nlev-20:nlev)
! if (lprnt) write(0,*)' icecldf=',icecldf(1,nlev-20:nlev)
! if (lprnt) write(0,*)' icldm=',icldm(1,nlev-20:nlev)
! if (lprnt) write(0,*)' qsfm=',qsfm(1,nlev-20:nlev)

  ! Initialize local variables

  ! local physical properties
  do k=1,nlev
    do i=1,mgncol
!     rho(i,k)    = p(i,k) / (r*t(i,k)*(one+fv*q(i,k)))
      rho(i,k)    = p(i,k) / (r*t(i,k))
      rhoinv(i,k) = one / rho(i,k)
      dv(i,k)     = 8.794E-5_r8 * t(i,k)**1.81_r8 / p(i,k)
      mu(i,k)     = 1.496E-6_r8 * t(i,k)*sqrt(t(i,k)) / (t(i,k) + 120._r8)
      sc(i,k)     = mu(i,k) / (rho(i,k)*dv(i,k))

  ! air density adjustment for fallspeed parameters
  ! includes air density correction factor to the
  ! power of 0.54 following Heymsfield and Bansemer 2007

      rhof(i,k) = (rhosu*rhoinv(i,k))**0.54_r8

      arn(i,k)  = ar*rhof(i,k)
      asn(i,k)  = as*rhof(i,k)
      acn(i,k)  = g*rhow/(18._r8*mu(i,k))
      tx1       = (rhosu*rhoinv(i,k))**0.35_r8
      ain(i,k)  = ai*tx1
      ajn(i,k)  = aj*tx1

  !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  ! Get humidity and saturation vapor pressures

! do k=1,nlev
!    do i=1,mgncol
!       relvar(i,k)      = relvar_i
        accre_enhan(i,k) = accre_enhan_i
!       call qsat_water(t(i,k), p(i,k), esl(i,k), qvl(i,k))
        esl(i,k)         = min(fpvsl(t(i,k)), p(i,k))
        qvl(i,k)         = epsqs*esl(i,k) / (p(i,k)-omeps*esl(i,k))


        ! make sure when above freezing that esi=esl, not active yet
        if (t(i,k) >= tmelt) then
           esi(i,k) = esl(i,k)
           qvi(i,k) = qvl(i,k)
        else
!          call qsat_ice(t(i,k), p(i,k), esi(i,k), qvi(i,k))
           esi(i,k) = min(fpvsi(t(i,k)), p(i,k))
           qvi(i,k) = epsqs*esi(i,k) / (p(i,k)-omeps*esi(i,k))
        end if

           ! Scale the water saturation values to reflect subgrid scale
           ! ice cloud fraction, where ice clouds begin forming at a
           ! gridbox average relative humidity of rhmini (not 1).
           !
           ! NOTE: For subcolumns and other non-subgrid clouds, qsfm will be 1.
        qvi(i,k) = qsfm(i,k) * qvi(i,k)
!       esi(i,k) = qsfm(i,k) * esi(i,k)
        qvl(i,k) = qsfm(i,k) * qvl(i,k)
!       esl(i,k) = qsfm(i,k) * esl(i,k)

        relhum(i,k) = max(zero, min(q(i,k)/max(qvl(i,k), qsmall), two))
     end do
  end do


  !===============================================

  ! set mtime here to avoid answer-changing
  mtime = deltat

  ! initialize microphysics output
  do k=1,nlev
    do i=1,mgncol
      qcsevap(i,k)     = zero
      qisevap(i,k)     = zero
      qvres(i,k)       = zero
      cmeitot(i,k)     = zero
      vtrmc(i,k)       = zero
      vtrmi(i,k)       = zero
      qcsedten(i,k)    = zero
      qisedten(i,k)    = zero
      qrsedten(i,k)    = zero
      qssedten(i,k)    = zero

      pratot(i,k)      = zero
      prctot(i,k)      = zero
      mnuccctot(i,k)   = zero
      mnuccttot(i,k)   = zero
      msacwitot(i,k)   = zero
      psacwstot(i,k)   = zero
      bergstot(i,k)    = zero
      bergtot(i,k)     = zero
      melttot(i,k)     = zero
      homotot(i,k)     = zero
      qcrestot(i,k)    = zero
      prcitot(i,k)     = zero
      praitot(i,k)     = zero
      qirestot(i,k)    = zero
      mnuccrtot(i,k)   = zero
      pracstot(i,k)    = zero
      meltsdttot(i,k)  = zero
      frzrdttot(i,k)   = zero
      mnuccdtot(i,k)   = zero

      rflx(i,k+1)      = zero
      sflx(i,k+1)      = zero
      lflx(i,k+1)      = zero
      iflx(i,k+1)      = zero

  ! initialize precip output

      qrout(i,k)       = zero
      qsout(i,k)       = zero
      nrout(i,k)       = zero
      nsout(i,k)       = zero

  ! for refl calc
      rainrt(i,k)      = zero

  ! initialize rain size
      rercld(i,k)      = zero

      qcsinksum_rate1ord(i,k) = zero

  ! initialize variables for trop_mozart
      nevapr(i,k)      = zero
      prer_evap(i,k)   = zero
      evapsnow(i,k)    = zero
      am_evp_st(i,k)   = zero
      prain(i,k)       = zero
      prodsnow(i,k)    = zero
      cmeout(i,k)      = zero

      precip_frac(i,k) = mincld

      lamc(i,k)        = zero

  ! initialize microphysical tendencies

      tlat(i,k)   = zero
      qvlat(i,k)  = zero
      qctend(i,k) = zero
      qitend(i,k) = zero
      qstend(i,k) = zero
      qrtend(i,k) = zero
      nctend(i,k) = zero
      nitend(i,k) = zero
      nrtend(i,k) = zero
      nstend(i,k) = zero

  ! initialize in-cloud and in-precip quantities to zero
      qcic(i,k)   = zero
      qiic(i,k)   = zero
      qsic(i,k)   = zero
      qric(i,k)   = zero

      ncic(i,k)   = zero
      niic(i,k)   = zero
      nsic(i,k)   = zero
      nric(i,k)   = zero

  ! initialize precip fallspeeds to zero
      ums(i,k)    = zero 
      uns(i,k)    = zero
      umr(i,k)    = zero
      unr(i,k)    = zero

  ! initialize limiter for output
      qcrat(i,k)  = one

  ! Many outputs have to be initialized here at the top to work around
  ! ifort problems, even if they are always overwritten later.
      effc(i,k)    = ten
      lamcrad(i,k) = zero
      pgamrad(i,k) = zero
      effc_fn(i,k) = ten
      effi(i,k)    = 25._r8
      sadice(i,k)  = zero
      sadsnow(i,k) = zero
      deffi(i,k)   = 50._r8

      qrout2(i,k) = zero
      nrout2(i,k) = zero
      drout2(i,k) = zero
      qsout2(i,k) = zero
      nsout2(i,k) = zero
      dsout(i,k)  = zero
      dsout2(i,k) = zero

      freqr(i,k) = zero
      freqs(i,k) = zero

      reff_rain(i,k) = zero
      reff_snow(i,k) = zero

      refl(i,k)   = -9999._r8
      arefl(i,k)  = zero
      areflz(i,k) = zero
      frefl(i,k)  = zero
      csrfl(i,k)  = zero
      acsrfl(i,k) = zero
      fcsrfl(i,k) = zero

      ncal(i,k)   = zero
      ncai(i,k)   = zero

      nfice(i,k)  = zero
      npccn(i,k)  = zero
    enddo
  enddo

  if(iccn == 1) then
    do k=1,nlev
      do i=1,mgncol
        npccn(i,k) = npccnin(i,k)
      enddo
    enddo
  else
    do k=1,nlev
      do i=1,mgncol
        npccn(i,k) = max((npccnin(i,k)*lcldm(i,k)-nc(i,k))*oneodt, zero)
      enddo
    enddo
  endif
  ! initialize precip at surface

  do i=1,mgncol
      prect(i) = zero
      preci(i) = zero
  enddo

  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
  ! droplet activation
  ! get provisional droplet number after activation. This is used for
  ! all microphysical process calculations, for consistency with update of
  ! droplet mass before microphysics

  ! calculate potential for droplet activation if cloud water is present
  ! tendency from activation (npccn) is read in from companion routine

  ! output activated liquid and ice (convert from #/kg -> #/m3)
  !--------------------------------------------------
  where (qc >= qsmall)
     npccn = max((npccnin*lcldm-nc)*oneodt, zero)
     nc    = max(nc + npccn*deltat, zero)
     ncal  = nc*rho/lcldm ! sghan minimum in #/cm3
  elsewhere
     ncal = zero
  end where

  if (iccn == 1) then
    do k=1,nlev
      do i=1,mgncol
        if (t(i,k) < icenuct) then
          ncai(i,k) = 0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k))) * 1000._r8
!         ncai(i,k) = min(ncai(i,k), 208.9e3_r8)
          ncai(i,k) = min(ncai(i,k), 355.0e3_r8)
          naai(i,k) = (ncai(i,k)*rhoinv(i,k) + naai(i,k))  * half
          ncai(i,k) = naai(i,k)*rho(i,k)
        else
          naai(i,k) = zero
          ncai(i,k) = zero
        endif
      enddo
    enddo
  elseif (iccn == 2) then
    do k=1,nlev
      do i=1,mgncol
        if (t(i,k) < icenuct) then
          ncai(i,k) = naai(i,k)*rho(i,k)
        else
          naai(i,k) = zero
          ncai(i,k) = zero
        endif
      enddo
    enddo
  else
    do k=1,nlev
      do i=1,mgncol
        if (t(i,k) < icenuct) then
          ncai(i,k) = 0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k))) * 1000._r8
          ncai(i,k) = min(ncai(i,k), 355.0e3_r8)
          naai(i,k) = ncai(i,k)*rhoinv(i,k)
        else
          naai(i,k) = zero
          ncai(i,k) = zero
        endif
      enddo
    enddo
    do k=1,nlev
      do i=1,mgncol
        naai(i,k) = zero
        ncai(i,k) = zero
      enddo
    enddo
  endif


  !===============================================

  ! ice nucleation if activated nuclei exist at t<-5C AND rhmini + 5%
  !
  ! NOTE: If using gridbox average values, condensation will not occur until rh=1,
  ! so the threshold seems like it should be 1.05 and not rhmini + 0.05. For subgrid
  ! clouds (using rhmini and qsfacm), the relhum has already been adjusted, and thus
  ! the nucleation threshold should also be 1.05 and not rhmini + 0.05.

  !-------------------------------------------------------

  if (do_cldice) then
     where (naai > zero .and. t < icenuct .and. relhum*esl/esi > 1.05_r8)

        !if NAAI > 0. then set numice = naai (as before)
        !note: this is gridbox averaged
        nnuccd = (naai-ni/icldm)/mtime*icldm
        nnuccd = max(nnuccd, zero)
        nimax  = naai*icldm

        !Calc mass of new particles using new crystal mass...
        !also this will be multiplied by mtime as nnuccd is...

        mnuccd = nnuccd * mi0

     elsewhere
        nnuccd = zero
        nimax  = zero
        mnuccd = zero
     end where

  end if


  !=============================================================================
  do k=1,nlev

     do i=1,mgncol

        ! calculate instantaneous precip processes (melting and homogeneous freezing)

        ! melting of snow at +2 C

        if (t(i,k) > snowmelt) then
           if (qs(i,k) > zero) then

              ! make sure melting snow doesn't reduce temperature below threshold
              dum = -(xlf/cpp) * qs(i,k)
              if (t(i,k)+dum < snowmelt) then
                 dum = min(one, max(zero, (cpp/xlf)*(t(i,k)-snowmelt)/qs(i,k)))
              else
                 dum = one
              end if

              minstsm(i,k) = dum*qs(i,k)
              ninstsm(i,k) = dum*ns(i,k)

              dum1 = - minstsm(i,k) * (xlf*oneodt)
              tlat(i,k)       = tlat(i,k)       + dum1
              meltsdttot(i,k) = meltsdttot(i,k) + dum1

              qs(i,k) = max(qs(i,k) - minstsm(i,k), zero)
              ns(i,k) = max(ns(i,k) - ninstsm(i,k), zero)
              qr(i,k) = max(qr(i,k) + minstsm(i,k), zero)
              nr(i,k) = max(nr(i,k) + ninstsm(i,k), zero)
           end if
        end if

     end do
  end do 
! if (lprnt) write(0,*)' tlat1=',tlat(1,:)*deltat

  do k=1,nlev
    do i=1,mgncol
        ! freezing of rain at -5 C

        if (t(i,k) < rainfrze) then

           if (qr(i,k) > zero) then

              ! make sure freezing rain doesn't increase temperature above threshold
              dum = (xlf/cpp) * qr(i,k)
              if (t(i,k)+dum > rainfrze) then
                 dum = -(t(i,k)-rainfrze) * (cpp/xlf)
                 dum = min(one, max(zero, dum/qr(i,k)))
              else
                 dum = one
              end if

              minstrf(i,k) = dum*qr(i,k)
              ninstrf(i,k) = dum*nr(i,k)

              ! heating tendency
              dum1 =  minstrf(i,k) * (xlf*oneodt)
              tlat(i,k)      = tlat(i,k)      + dum1
              frzrdttot(i,k) = frzrdttot(i,k) + dum1

              qr(i,k) = max(qr(i,k) - minstrf(i,k), zero)
              nr(i,k) = max(nr(i,k) - ninstrf(i,k), zero)
              qs(i,k) = max(qs(i,k) + minstrf(i,k), zero)
              ns(i,k) = max(ns(i,k) + ninstrf(i,k), zero)

           end if
        end if
     end do
  end do 

! if (lprnt) write(0,*)' tlat2=',tlat(1,:)*deltat
  do k=1,nlev
    do i=1,mgncol
        ! obtain in-cloud values of cloud water/ice mixing ratios and number concentrations
        !-------------------------------------------------------
        ! for microphysical process calculations
        ! units are kg/kg for mixing ratio, 1/kg for number conc

        if (qc(i,k) >= qsmall) then
           dum = one / lcldm(i,k)
!          qcic(i,k) = min(qc(i,k)*dum, 5.e-3_r8) ! limit in-cloud values to 0.005 kg/kg
           qcic(i,k) = min(qc(i,k)*dum, 0.05_r8)  ! limit in-cloud values to 0.05 kg/kg
           ncic(i,k) = max(nc(i,k)*dum, zero)

           ! specify droplet concentration
           if (nccons) then
              ncic(i,k) = ncnst * rhoinv(i,k)
           end if
        else
           qcic(i,k) = zero
           ncic(i,k) = zero
        end if

        if (qi(i,k) >= qsmall) then
           dum = one / icldm(i,k)
!          qiic(i,k) = min(qi(i,k)*dum, 5.e-3_r8) ! limit in-cloud values to 0.005 kg/kg
           qiic(i,k) = min(qi(i,k)*dum, 0.05_r8)  ! limit in-cloud values to 0.05 kg/kg
           niic(i,k) = max(ni(i,k)*dum, zero)

           ! switch for specification of cloud ice number
           if (nicons) then
              niic(i,k) = ninst * rhoinv(i,k)
           end if
        else
           qiic(i,k) = zero
           niic(i,k) = zero
        end if

     end do
  end do

  !========================================================================

  ! for sub-columns cldm has already been set to 1 if cloud
  ! water or ice is present, so precip_frac will be correctly set below
  ! and nothing extra needs to be done here

  precip_frac = cldm

  micro_vert_loop: do k=1,nlev

     if (trim(micro_mg_precip_frac_method) == 'in_cloud') then

        if (k /= 1) then
           where (qc(:,k) < qsmall .and. qi(:,k) < qsmall)
              precip_frac(:,k) = precip_frac(:,k-1)
           end where
        endif

     else if (trim(micro_mg_precip_frac_method) == 'max_overlap') then

        ! calculate precip fraction based on maximum overlap assumption

        ! if rain or snow mix ratios are smaller than threshold,
        ! then leave precip_frac as cloud fraction at current level
        if (k /= 1) then
           where (qr(:,k-1) >= qsmall .or. qs(:,k-1) >= qsmall)
              precip_frac(:,k) = max(precip_frac(:,k-1),precip_frac(:,k))
           end where
        end if

     endif


     !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
     ! get size distribution parameters based on in-cloud cloud water
     ! these calculations also ensure consistency between number and mixing ratio
     !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

     ! cloud liquid
     !-------------------------------------------

     call size_dist_param_liq(mg_liq_props, qcic(:,k), ncic(:,k), rho(:,k), &
                              pgam(:,k), lamc(:,k), mgncol)


     !========================================================================
     ! autoconversion of cloud liquid water to rain
     ! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc
     ! minimum qc of 1 x 10^-8 prevents floating point error

     if (.not. do_sb_physics) then
        call kk2000_liq_autoconversion(microp_uniform, qcic(:,k), &
             ncic(:,k), rho(:,k), relvar(:,k), prc(:,k), nprc(:,k), nprc1(:,k), mgncol)
     endif

     ! assign qric based on prognostic qr, using assumed precip fraction
     ! note: this could be moved above for consistency with qcic and qiic calculations
     do i=1,mgncol
       if (precip_frac(i,k) > mincld) then
         dum = one / precip_frac(i,k)
       else
         dum = zero
       endif
!      qric(i,k) = min(qr(i,k)*dum, 0.01_r8) ! limit in-precip mixing ratios to 10 g/kg
       qric(i,k) = min(qr(i,k)*dum, 0.05_r8) ! limit in-precip mixing ratios to 50 g/kg
       nric(i,k) = nr(i,k) * dum


     ! add autoconversion to precip from above to get provisional rain mixing ratio
     ! and number concentration (qric and nric)

       if(qric(i,k) < qsmall) then
         qric(i,k) = zero
         nric(i,k) = zero
       endif

     ! make sure number concentration is a positive number to avoid
     ! taking root of negative later

       nric(i,k) = max(nric(i,k),zero)
     enddo 
     ! Get size distribution parameters for cloud ice

     call size_dist_param_ice(mg_ice_props, qiic(:,k), niic(:,k), &
                              lami(:,k), mgncol, n0=n0i(:,k))

!    call size_dist_param_basic(mg_ice_props, qiic(:,k), niic(:,k), &
!         lami(:,k), mgncol, n0=n0i(:,k))
  
     ! Alternative autoconversion 
     if (do_sb_physics) then
       if  (do_liq_liu) then
         call liu_liq_autoconversion(pgam(:,k),qcic(:,k),ncic(:,k), &
              qric(:,k),rho(:,k),relvar(:,k),prc(:,k),nprc(:,k),nprc1(:,k),mgncol)
       else
         call sb2001v2_liq_autoconversion(pgam(:,k),qcic(:,k),ncic(:,k), &
              qric(:,k),rho(:,k),relvar(:,k),prc(:,k),nprc(:,k),nprc1(:,k), mgncol)
       endif
     endif

     !.......................................................................
     ! Autoconversion of cloud ice to snow
     ! similar to Ferrier (1994)

     if (do_cldice) then
        do i=1,mgncol
          if (qiic(i,k) >= qimax) then
!         if (qi(i,k) >= qimax) then
            ts_au_loc(i) = ts_au_min
          elseif (qiic(i,k) <= qimin) then
!         elseif (qi(i,k) <= qimin) then
            ts_au_loc(i) = ts_au
          else
!           ts_au_loc(i) = (ts_au*(qimax-qi(i,k)) + ts_au_min*(qi(i,k)-qimin)) * qiinv
            ts_au_loc(i) = (ts_au*(qimax-qiic(i,k)) + ts_au_min*(qiic(i,k)-qimin)) * qiinv
!           ts_au_loc(i) = ts_au * exp(-tsfac*(qiic(i,k)-qimin))
          endif
        enddo

        if(do_ice_gmao) then
          call gmao_ice_autoconversion(t(:,k), qiic(:,k), niic(:,k), lami(:,k),     &
               n0i(:,k), dcs, ts_au_loc(:), prci(:,k), nprci(:,k), mgncol)
        else
          call ice_autoconversion(t(:,k), qiic(:,k), lami(:,k), n0i(:,k),           &
                         dcs, ts_au_loc(:), prci(:,k), nprci(:,k), mgncol)
        end if
     !else
        ! Add in the particles that we have already converted to snow, and
        ! don't do any further autoconversion of ice.
        !prci(:,k)  = tnd_qsnow(:,k) / cldm(:,k)
        !nprci(:,k) = tnd_nsnow(:,k) / cldm(:,k)
     end if

     ! note, currently we don't have this
     ! inside the do_cldice block, should be changed later
     ! assign qsic based on prognostic qs, using assumed precip fraction
     do i=1,mgncol
       if (precip_frac(i,k) > mincld) then
         dum = one / precip_frac(i,k)
       else
         dum = zero
       endif
!      qsic(i,k) = min(qs(i,k)*dum, 0.01_r8) ! limit in-precip mixing ratios to 10 g/kg
!      qsic(i,k) = min(qs(i,k)*dum, 0.05_r8) ! limit in-precip mixing ratios to 50 g/kg
       qsic(i,k) = min(qs(i,k)*dum, 0.10_r8) ! limit in-precip mixing ratios to 50 g/kg
       nsic(i,k) = ns(i,k) * dum

     ! if precip mix ratio is zero so should number concentration

       if(qsic(i,k) < qsmall) then
         qsic(i,k) = zero
         nsic(i,k) = zero
       endif

     ! make sure number concentration is a positive number to avoid
     ! taking root of negative later

       nsic(i,k) = max(nsic(i,k), zero)
     enddo

     !.......................................................................
     ! get size distribution parameters for precip
     !......................................................................
     ! rain

     call size_dist_param_basic(mg_rain_props, qric(:,k), nric(:,k), &
          lamr(:,k), mgncol, n0=n0r(:,k))

     do i=1,mgncol
       if (lamr(i,k) >= qsmall) then
          dum  = arn(i,k) / lamr(i,k)**br
          dum1 = 9.1_r8*rhof(i,k)

        ! provisional rain number and mass weighted mean fallspeed (m/s)

         umr(i,k) = min(dum1, dum*gamma_br_plus4*oneo6)
         unr(i,k) = min(dum1, dum*gamma_br_plus1)
       else

         umr(i,k) = zero
         unr(i,k) = zero
       endif
     enddo

     !......................................................................
     ! snow

     call size_dist_param_basic(mg_snow_props, qsic(:,k), nsic(:,k), &
                                lams(:,k), mgncol, n0=n0s(:,k))

     do i=1,mgncol
       if (lams(i,k) >= qsmall) then

        ! provisional snow number and mass weighted mean fallspeed (m/s)

         dum  = asn(i,k) / lams(i,k)**bs
         dum1 = 1.2_r8*rhof(i,k)
         ums(i,k) = min(dum1, dum*gamma_bs_plus4*oneo6)
         uns(i,k) = min(dum1, dum*gamma_bs_plus1)

       else
         ums(i,k) = zero
         uns(i,k) = zero
       endif
     enddo

     if (do_cldice) then
        if (.not. use_hetfrz_classnuc) then

           ! heterogeneous freezing of cloud water
           !----------------------------------------------

           call immersion_freezing(microp_uniform, t(:,k), pgam(:,k), lamc(:,k), &
                qcic(:,k), ncic(:,k), relvar(:,k), mnuccc(:,k), nnuccc(:,k), mgncol)

           ! make sure number of droplets frozen does not exceed available ice nuclei concentration
           ! this prevents 'runaway' droplet freezing

           where (qcic(:,k) >= qsmall .and. t(:,k) < 269.15_r8)
              where (nnuccc(:,k)*lcldm(:,k) > nnuccd(:,k))
                 ! scale mixing ratio of droplet freezing with limit
                 mnuccc(:,k) = mnuccc(:,k)*(nnuccd(:,k)/(nnuccc(:,k)*lcldm(:,k)))
                 nnuccc(:,k) = nnuccd(:,k)/lcldm(:,k)
              end where
           end where

           mdust = size(rndst,3)
           call contact_freezing(microp_uniform, t(:,k), p(:,k), rndst(:,k,:), &
                nacon(:,k,:), pgam(:,k), lamc(:,k), qcic(:,k), ncic(:,k), &
                relvar(:,k), mnucct(:,k), nnucct(:,k), mgncol, mdust)

           mnudep(:,k) = zero
           nnudep(:,k) = zero

        !else

           ! Mass of droplets frozen is the average droplet mass, except
           ! with two limiters: concentration must be at least 1/cm^3, and
           ! mass must be at least the minimum defined above.
           !mi0l = qcic(:,k)/max(ncic(:,k), 1.0e6_r8/rho(:,k))
           !mi0l = max(mi0l_min, mi0l)

           !where (qcic(:,k) >= qsmall)
              !nnuccc(:,k) = frzimm(:,k)*1.0e6_r8/rho(:,k)
              !mnuccc(:,k) = nnuccc(:,k)*mi0l

              !nnucct(:,k) = frzcnt(:,k)*1.0e6_r8/rho(:,k)
              !mnucct(:,k) = nnucct(:,k)*mi0l

              !nnudep(:,k) = frzdep(:,k)*1.0e6_r8/rho(:,k)
              !mnudep(:,k) = nnudep(:,k)*mi0
           !elsewhere
              !nnuccc(:,k) = 0._r8
              !mnuccc(:,k) = 0._r8

              !nnucct(:,k) = 0._r8
              !mnucct(:,k) = 0._r8

              !nnudep(:,k) = 0._r8
              !mnudep(:,k) = 0._r8
           !end where

        end if

     else
        do i=1,mgncol
          mnuccc(i,k) = zero
          nnuccc(i,k) = zero
          mnucct(i,k) = zero
          nnucct(i,k) = zero
          mnudep(i,k) = zero
          nnudep(i,k) = zero
        enddo
     end if

     call snow_self_aggregation(t(:,k), rho(:,k), asn(:,k), rhosn, qsic(:,k), nsic(:,k), &
          nsagg(:,k), mgncol)

     call accrete_cloud_water_snow(t(:,k), rho(:,k), asn(:,k), uns(:,k), mu(:,k), &
          qcic(:,k), ncic(:,k), qsic(:,k), pgam(:,k), lamc(:,k), lams(:,k), n0s(:,k), &
          psacws(:,k), npsacws(:,k), mgncol)

     if (do_cldice) then
        call secondary_ice_production(t(:,k), psacws(:,k), msacwi(:,k), nsacwi(:,k), mgncol)
     else
        nsacwi(:,k) = zero
        msacwi(:,k) = zero
     end if

     call accrete_rain_snow(t(:,k), rho(:,k), umr(:,k), ums(:,k), unr(:,k), uns(:,k), &
          qric(:,k), qsic(:,k), lamr(:,k), n0r(:,k), lams(:,k), n0s(:,k), &
          pracs(:,k), npracs(:,k), mgncol)

     call heterogeneous_rain_freezing(t(:,k), qric(:,k), nric(:,k), lamr(:,k), &
          mnuccr(:,k), nnuccr(:,k), mgncol)

     if (do_sb_physics) then
       call sb2001v2_accre_cld_water_rain(qcic(:,k), ncic(:,k), qric(:,k), &
            rho(:,k), relvar(:,k), pra(:,k), npra(:,k), mgncol)     
     else
       call accrete_cloud_water_rain(microp_uniform, qric(:,k), qcic(:,k), &
            ncic(:,k), relvar(:,k), accre_enhan(:,k), pra(:,k), npra(:,k), mgncol)
     endif

     call self_collection_rain(rho(:,k), qric(:,k), nric(:,k), nragg(:,k), mgncol)

     if (do_cldice) then
        call accrete_cloud_ice_snow(t(:,k), rho(:,k), asn(:,k), qiic(:,k), niic(:,k), &
             qsic(:,k), lams(:,k), n0s(:,k), prai(:,k), nprai(:,k), mgncol)
     else
        prai(:,k)  = zero
        nprai(:,k) = zero
     end if

     call evaporate_sublimate_precip(t(:,k), rho(:,k), &
          dv(:,k), mu(:,k), sc(:,k), q(:,k), qvl(:,k), qvi(:,k), &
          lcldm(:,k), precip_frac(:,k), arn(:,k), asn(:,k), qcic(:,k), qiic(:,k), &
          qric(:,k), qsic(:,k), lamr(:,k), n0r(:,k), lams(:,k), n0s(:,k), &
          pre(:,k), prds(:,k), am_evp_st(:,k), mgncol)

     call bergeron_process_snow(t(:,k), rho(:,k), dv(:,k), mu(:,k), sc(:,k), &
          qvl(:,k), qvi(:,k), asn(:,k), qcic(:,k), qsic(:,k), lams(:,k), n0s(:,k), &
          bergs(:,k), mgncol)

     bergs(:,k)=bergs(:,k)*micro_mg_berg_eff_factor

     !+++PMC 12/3/12 - NEW VAPOR DEP/SUBLIMATION GOES HERE!!!
     if (do_cldice) then

       call ice_deposition_sublimation(t(:,k), q(:,k), qi(:,k), ni(:,k), &
            icldm(:,k), rho(:,k), dv(:,k), qvl(:,k), qvi(:,k), &
            berg(:,k), vap_dep(:,k), ice_sublim(:,k), mgncol)

       do i=1,mgncol
! sublimation should not exceed available ice
         ice_sublim(i,k) = max(ice_sublim(i,k), -qi(i,k)*oneodt)

         berg(i,k) = berg(i,k)*micro_mg_berg_eff_factor

         if (vap_dep(i,k) < zero .and. qi(i,k) > qsmall .and. icldm(i,k) > mincld) then
           nsubi(i,k) = vap_dep(i,k) * ni(i,k) / (qi(i,k) * icldm(i,k))
         else
           nsubi(i,k) = zero
         endif

        ! bergeron process should not reduce nc unless
        ! all ql is removed (which is handled elsewhere)
        !in fact, nothing in this entire file makes nsubc nonzero.

         nsubc(i,k) = zero
       enddo

     end if !do_cldice
     !---PMC 12/3/12

     do i=1,mgncol

        ! conservation to ensure no negative values of cloud water/precipitation
        ! in case microphysical process rates are large
        !===================================================================

        ! note: for check on conservation, processes are multiplied by omsm
        ! to prevent problems due to round off error

        ! conservation of qc
        !-------------------------------------------------------------------

        dum = ((prc(i,k)+pra(i,k)+mnuccc(i,k)+mnucct(i,k)+msacwi(i,k)+ &
                psacws(i,k)+bergs(i,k))*lcldm(i,k)+berg(i,k))*deltat

        if (dum > qc(i,k)) then
           ratio = qc(i,k) / dum * omsm

           prc(i,k)    = ratio * prc(i,k)
           pra(i,k)    = ratio * pra(i,k)
           mnuccc(i,k) = ratio * mnuccc(i,k)
           mnucct(i,k) = ratio * mnucct(i,k)
           msacwi(i,k) = ratio * msacwi(i,k)
           psacws(i,k) = ratio * psacws(i,k)
           bergs(i,k)  = ratio * bergs(i,k)
           berg(i,k)   = ratio * berg(i,k)
           qcrat(i,k)  = ratio
        else
           qcrat(i,k) = one
        end if

        !PMC 12/3/12: ratio is also frac of step w/ liquid.
        !thus we apply berg for "ratio" of timestep and vapor
        !deposition for the remaining frac of the timestep.
        if (qc(i,k) >= qsmall) then
           vap_dep(i,k) = vap_dep(i,k)*(1._r8-qcrat(i,k))
        end if

     end do

     do i=1,mgncol

        !=================================================================
        ! apply limiter to ensure that ice/snow sublimation and rain evap
        ! don't push conditions into supersaturation, and ice deposition/nucleation don't
        ! push conditions into sub-saturation
        ! note this is done after qc conservation since we don't know how large
        ! vap_dep is before then
        ! estimates are only approximate since other process terms haven't been limited
        ! for conservation yet

        ! first limit ice deposition/nucleation vap_dep + mnuccd

        dum1 = vap_dep(i,k) + mnuccd(i,k)
        if (dum1 > 1.e-20_r8) then
           dum = (q(i,k)-qvi(i,k))/(one + xxls_squared*qvi(i,k)/(cpp*rv*t(i,k)*t(i,k)))*oneodt
           dum = max(dum, zero)
           if (dum1 > dum) then
              ! Allocate the limited "dum" tendency to mnuccd and vap_dep
              ! processes. Don't divide by cloud fraction; these are grid-
              ! mean rates.
              dum1 = mnuccd(i,k) / (vap_dep(i,k)+mnuccd(i,k))
              mnuccd(i,k)  = dum*dum1
              vap_dep(i,k) = dum - mnuccd(i,k)
           end if
        end if

     end do

     do i=1,mgncol

        !===================================================================
        ! conservation of nc
        !-------------------------------------------------------------------
        dum = (nprc1(i,k)+npra(i,k)+nnuccc(i,k)+nnucct(i,k)+ &
               npsacws(i,k)-nsubc(i,k))*lcldm(i,k) * deltat

        if (dum > nc(i,k)) then
           ratio = nc(i,k) / dum * omsm

           nprc1(i,k)   = ratio * nprc1(i,k)
           npra(i,k)    = ratio * npra(i,k)
           nnuccc(i,k)  = ratio * nnuccc(i,k)
           nnucct(i,k)  = ratio * nnucct(i,k)
           npsacws(i,k) = ratio * npsacws(i,k)
           nsubc(i,k)   = ratio * nsubc(i,k)
        endif

        mnuccri(i,k) = zero
        nnuccri(i,k) = zero

        if (do_cldice) then

           ! freezing of rain to produce ice if mean rain size is smaller than Dcs
           if (lamr(i,k) > qsmall) then
             if (one/lamr(i,k) < Dcs) then
               mnuccri(i,k) = mnuccr(i,k)
               nnuccri(i,k) = nnuccr(i,k)
               mnuccr(i,k)  = zero 
               nnuccr(i,k)  = zero
             endif
           endif
        endif

     enddo

     do i=1,mgncol

        ! conservation of rain mixing ratio
        !-------------------------------------------------------------------
        dum1 = -pre(i,k) + pracs(i,k) + mnuccr(i,k) + mnuccri(i,k)
        dum3 = dum1 * precip_frac(i,k)
        dum2 = (pra(i,k)+prc(i,k))*lcldm(i,k)
        dum = (dum3 - dum2) * deltat

        ! note that qrtend is included below because of instantaneous freezing/melt
        if (dum > qr(i,k) .and. dum1 >= qsmall) then
           ratio = (qr(i,k)*oneodt + dum2) / dum3 * omsm
           pre(i,k)     = ratio * pre(i,k)
           pracs(i,k)   = ratio * pracs(i,k)
           mnuccr(i,k)  = ratio * mnuccr(i,k)
           mnuccri(i,k) = ratio * mnuccri(i,k)
        end if

     end do

     do i=1,mgncol

        ! conservation of rain number
        !-------------------------------------------------------------------

        ! Add evaporation of rain number.
        if (pre(i,k) < zero) then
           dum        = max(-one, pre(i,k)*deltat/qr(i,k))
           nsubr(i,k) = dum*nr(i,k) * oneodt
        else
           nsubr(i,k) = zero
        end if

     end do

     do i=1,mgncol

        dum1 = (-nsubr(i,k)+npracs(i,k)+nnuccr(i,k)+nnuccri(i,k)-nragg(i,k))*precip_frac(i,k)
        dum2 = nprc(i,k)*lcldm(i,k)
        dum = (dum1 - dum2) * deltat

        if (dum > nr(i,k)) then
           ratio = (nr(i,k)*oneodt + dum2) / dum1 * omsm

           nragg(i,k)   = ratio * nragg(i,k)
           npracs(i,k)  = ratio * npracs(i,k)
           nnuccr(i,k)  = ratio * nnuccr(i,k)
           nsubr(i,k)   = ratio * nsubr(i,k)
           nnuccri(i,k) = ratio * nnuccri(i,k)
        end if

     end do

     if (do_cldice) then

        do i=1,mgncol

           ! conservation of qi
           !-------------------------------------------------------------------

           dum1 = (prci(i,k)+prai(i,k))*icldm(i,k)-ice_sublim(i,k)
           dum2 = vap_dep(i,k)+berg(i,k)+mnuccd(i,k)                       &
                + (mnuccc(i,k)+mnucct(i,k)+mnudep(i,k)+msacwi(i,k))*lcldm(i,k) &
                + mnuccri(i,k)*precip_frac(i,k)
           dum  = (dum1 - dum2) * deltat

           if (dum > qi(i,k)) then
              ratio = (qi(i,k)*oneodt + dum2) / dum1 * omsm

              prci(i,k)       = ratio * prci(i,k)
              prai(i,k)       = ratio * prai(i,k)
              ice_sublim(i,k) = ratio * ice_sublim(i,k)
           end if

        end do

     end if

     if (do_cldice) then

        do i=1,mgncol

           ! conservation of ni
           !-------------------------------------------------------------------
           if (use_hetfrz_classnuc) then
              tmpfrz = nnuccc(i,k)
           else
              tmpfrz = zero
           end if
           dum1 = (nprci(i,k)+nprai(i,k)-nsubi(i,k))*icldm(i,k)
           dum2 = nnuccd(i,k)+(nnucct(i,k)+tmpfrz+nnudep(i,k)+nsacwi(i,k))*lcldm(i,k) &
                             + nnuccri(i,k)*precip_frac(i,k)
           dum  = (dum1 - dum2) * deltat

           if (dum > ni(i,k)) then
              ratio = (ni(i,k)*oneodt + dum2) / dum1 * omsm

              nprci(i,k) = ratio * nprci(i,k)
              nprai(i,k) = ratio * nprai(i,k)
              nsubi(i,k) = ratio * nsubi(i,k)
           end if

        end do

     end if

     do i=1,mgncol

        ! conservation of snow mixing ratio
        !-------------------------------------------------------------------
        dum1 = - prds(i,k) * precip_frac(i,k)
        dum2 = (pracs(i,k)+mnuccr(i,k))*precip_frac(i,k)                            &
             + (prai(i,k)+prci(i,k))*icldm(i,k) + (bergs(i,k)+psacws(i,k))*lcldm(i,k)
        dum  = (dum1 - dum2) * deltat

        if (dum > qs(i,k) .and. -prds(i,k) >= qsmall) then
           ratio = (qs(i,k)*oneodt + dum2) / dum1 * omsm

           prds(i,k) = ratio * prds(i,k)
        end if

     end do

     do i=1,mgncol

        ! conservation of snow number
        !-------------------------------------------------------------------
        ! calculate loss of number due to sublimation
        ! for now neglect sublimation of ns
        nsubs(i,k) = zero

        dum1 = precip_frac(i,k)* (-nsubs(i,k)-nsagg(i,k))
        dum2 = nnuccr(i,k)*precip_frac(i,k) + nprci(i,k)*icldm(i,k)
        dum  = (dum1 - dum2) * deltat

        if (dum > ns(i,k)) then
           ratio = (ns(i,k)*oneodt + dum2) / dum1 * omsm

           nsubs(i,k) = ratio * nsubs(i,k)
           nsagg(i,k) = ratio * nsagg(i,k)
        end if

     end do

     do i=1,mgncol

        ! next limit ice and snow sublimation and rain evaporation
        ! get estimate of q and t at end of time step
        ! don't include other microphysical processes since they haven't
        ! been limited via conservation checks yet

        tx1 = pre(i,k)  * precip_frac(i,k)
        tx2 = prds(i,k) * precip_frac(i,k)
        tx3 = tx1 + tx2 + ice_sublim(i,k)
        if (tx3 < -1.e-20_r8) then

           tx4  = tx2 + ice_sublim(i,k) + vap_dep(i,k) + mnuccd(i,k)
           qtmp = q(i,k) - (tx1 + tx4) * deltat
           ttmp = t(i,k) + (tx1*xxlv + tx4*xxls) * (deltat/cpp)

           ! use rhw to allow ice supersaturation
           ! call qsat_water(ttmp, p(i,k), esn, qvn)
           esn = min(fpvsl(ttmp), p(i,k))
           qvn = epsqs*esn/(p(i,k)-omeps*esn) * qsfm(i,k)

           ! modify ice/precip evaporation rate if q > qsat
           if (qtmp > qvn) then

              tx4  = one / tx3
              dum1 = tx1 * tx4
              dum2 = tx2 * tx4
              ! recalculate q and t after vap_dep and mnuccd but without evap or sublim
              tx5  = (vap_dep(i,k)+mnuccd(i,k)) * deltat
              qtmp = q(i,k) - tx5
              ttmp = t(i,k) + tx5 * (xxls/cpp)

              ! use rhw to allow ice supersaturation
              !call qsat_water(ttmp, p(i,k), esn, qvn)
              esn = min(fpvsl(ttmp), p(i,k))
              qvn = epsqs*esn / (p(i,k)-omeps*esn) * qsfm(i,k)

              dum = (qtmp-qvn) / (one + xxlv_squared*qvn/(cpp*rv*ttmp*ttmp))
              dum = min(dum, zero)

              ! modify rates if needed, divide by precip_frac to get local (in-precip) value
              if (precip_frac(i,k) > mincld) then
                tx4 = oneodt / precip_frac(i,k)
              else
                tx4 = zero
              endif
              pre(i,k) = dum*dum1*tx4

              ! do separately using RHI for prds and ice_sublim
              !call qsat_ice(ttmp, p(i,k), esn, qvn)
              esn = min(fpvsi(ttmp), p(i,k))
              qvn = epsqs*esn / (p(i,k)-omeps*esn) * qsfm(i,k)


              dum = (qtmp-qvn) / (one + xxls_squared*qvn/(cpp*rv*ttmp*ttmp))
              dum = min(dum, zero)

              ! modify rates if needed, divide by precip_frac to get local (in-precip) value
              prds(i,k) = dum*dum2*tx4

              ! don't divide ice_sublim by cloud fraction since it is grid-averaged
              dum1 = one - dum1 - dum2
              ice_sublim(i,k) = dum*dum1*oneodt
           end if
        end if

     end do

     ! Big "administration" loop enforces conservation, updates variables
     ! that accumulate over substeps, and sets output variables.

     do i=1,mgncol

        ! get tendencies due to microphysical conversion processes
        !==========================================================
        ! note: tendencies are multiplied by appropriate cloud/precip
        ! fraction to get grid-scale values
        ! note: vap_dep is already grid-average values

        ! The net tendencies need to be added to rather than overwritten,
        ! because they may have a value already set for instantaneous
        ! melting/freezing.


        qvlat(i,k) = qvlat(i,k) - (pre(i,k)+prds(i,k))*precip_frac(i,k)-&
             vap_dep(i,k)-ice_sublim(i,k)-mnuccd(i,k)-mnudep(i,k)*lcldm(i,k)


!     if (lprnt .and. k >= 60 ) &
!               write(0,*)' k=',k,' tlat=',tlat(i,k),' pre=',pre(i,k),' precip_frac=',precip_frac(i,k),&
!       ' prds=',prds(i,k),' vap_dep=',vap_dep(i,k),' ice_sublim=',ice_sublim(i,k), &
!       ' mnuccd=',mnuccd(i,k),' mnudep=',mnudep(i,k),' lcldm=',lcldm(i,k),' bergs=',bergs(i,k),       &
!       ' psacws=',psacws(i,k),' mnuccc=',mnuccc(i,k),' mnucct=',mnucct(i,k),' msacwi=',msacwi(i,k),   &
!       ' mnuccr=',mnuccr(i,k),   &
!       ' pracs=',pracs(i,k),' mnuccri=',mnuccri(i,k),' xlf=',xlf,' xxlv=',xxlv,' xxls=',xxls

        tlat(i,k) = tlat(i,k) + ((pre(i,k)*precip_frac(i,k)) &
             *xxlv+(prds(i,k)*precip_frac(i,k)+vap_dep(i,k)+ice_sublim(i,k)+mnuccd(i,k)+mnudep(i,k)*lcldm(i,k))*xxls+ &
             ((bergs(i,k)+psacws(i,k)+mnuccc(i,k)+mnucct(i,k)+msacwi(i,k))*lcldm(i,k)+(mnuccr(i,k)+ &
               pracs(i,k)+mnuccri(i,k))*precip_frac(i,k)+berg(i,k))*xlf)

!     if (lprnt .and. k >= 60) write(0,*)' k=',k,' tlat=',tlat(i,k)

        qctend(i,k) = qctend(i,k) + (-pra(i,k)-prc(i,k)-mnuccc(i,k)-mnucct(i,k)-msacwi(i,k)- &
                                      psacws(i,k)-bergs(i,k))*lcldm(i,k)-berg(i,k)

        if (do_cldice) then
           qitend(i,k) = qitend(i,k) + &
                (mnuccc(i,k)+mnucct(i,k)+mnudep(i,k)+msacwi(i,k))*lcldm(i,k)+(-prci(i,k)- &
                 prai(i,k))*icldm(i,k)+vap_dep(i,k)+berg(i,k)+ice_sublim(i,k)+ &
                 mnuccd(i,k)+mnuccri(i,k)*precip_frac(i,k)
        end if

        qrtend(i,k) = qrtend(i,k) + (pra(i,k)+prc(i,k))*lcldm(i,k)+(pre(i,k)-pracs(i,k)- &
                                     mnuccr(i,k)-mnuccri(i,k))*precip_frac(i,k)

        qstend(i,k) = qstend(i,k) + (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k) &
                                  + (prds(i,k)+pracs(i,k)+mnuccr(i,k))*precip_frac(i,k)


        cmeout(i,k) = vap_dep(i,k) + ice_sublim(i,k) + mnuccd(i,k)

        ! add output for cmei (accumulate)
        cmeitot(i,k) = vap_dep(i,k) + ice_sublim(i,k) + mnuccd(i,k)

        ! assign variables for trop_mozart, these are grid-average
        !-------------------------------------------------------------------
        ! evaporation/sublimation is stored here as positive term

        evapsnow(i,k)  = -prds(i,k) * precip_frac(i,k)
        nevapr(i,k)    = -pre(i,k)  * precip_frac(i,k)
        prer_evap(i,k) = -pre(i,k)  * precip_frac(i,k)

        ! change to make sure prain is positive: do not remove snow from
        ! prain used for wet deposition
        prain(i,k) = (pra(i,k)+prc(i,k))*lcldm(i,k)+(-pracs(i,k)- &
             mnuccr(i,k)-mnuccri(i,k))*precip_frac(i,k)
        prodsnow(i,k) = (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k)+(&
             pracs(i,k)+mnuccr(i,k))*precip_frac(i,k)

        ! following are used to calculate 1st order conversion rate of cloud water
        !    to rain and snow (1/s), for later use in aerosol wet removal routine
        ! previously, wetdepa used (prain/qc) for this, and the qc in wetdepa may be smaller than the qc
        !    used to calculate pra, prc, ... in this routine
        ! qcsinksum_rate1ord = { rate of direct transfer of cloud water to rain & snow }
        !                      (no cloud ice or bergeron terms)
        qcsinksum_rate1ord(i,k) = (pra(i,k)+prc(i,k)+psacws(i,k))*lcldm(i,k)
        ! Avoid zero/near-zero division.
        qcsinksum_rate1ord(i,k) = qcsinksum_rate1ord(i,k) / max(qc(i,k),1.0e-30_r8)


        ! microphysics output, note this is grid-averaged
        pratot(i,k)    = pra(i,k)    * lcldm(i,k)
        prctot(i,k)    = prc(i,k)    * lcldm(i,k)
        mnuccctot(i,k) = mnuccc(i,k) * lcldm(i,k)
        mnuccttot(i,k) = mnucct(i,k) * lcldm(i,k)
        msacwitot(i,k) = msacwi(i,k) * lcldm(i,k)
        psacwstot(i,k) = psacws(i,k) * lcldm(i,k)
        bergstot(i,k)  = bergs(i,k)  * lcldm(i,k)
        bergtot(i,k)   = berg(i,k)
        prcitot(i,k)   = prci(i,k)   * icldm(i,k)
        praitot(i,k)   = prai(i,k)   * icldm(i,k)
        mnuccdtot(i,k) = mnuccd(i,k) * icldm(i,k)

        pracstot(i,k)  = pracs(i,k)  * precip_frac(i,k)
        mnuccrtot(i,k) = mnuccr(i,k) * precip_frac(i,k)


        nctend(i,k) = nctend(i,k) + (-nnuccc(i,k)-nnucct(i,k)-npsacws(i,k)+nsubc(i,k) &
                                    - npra(i,k)-nprc1(i,k))*lcldm(i,k)

        if (do_cldice) then
           if (use_hetfrz_classnuc) then
              tmpfrz = nnuccc(i,k)
           else
              tmpfrz = zero
           end if
           nitend(i,k) = nitend(i,k) + nnuccd(i,k)+ &
                (nnucct(i,k)+tmpfrz+nnudep(i,k)+nsacwi(i,k))*lcldm(i,k)+(nsubi(i,k)-nprci(i,k)- &
                 nprai(i,k))*icldm(i,k)+nnuccri(i,k)*precip_frac(i,k)
        end if

        nstend(i,k) = nstend(i,k) + (nsubs(i,k)+nsagg(i,k)+nnuccr(i,k))*precip_frac(i,k) &
                                  + nprci(i,k)*icldm(i,k)

        nrtend(i,k) = nrtend(i,k) + nprc(i,k)*lcldm(i,k)+(nsubr(i,k)-npracs(i,k)-nnuccr(i,k) &
                                  - nnuccri(i,k)+nragg(i,k))*precip_frac(i,k)

        ! make sure that ni at advanced time step does not exceed
        ! maximum (existing N + source terms*dt), which is possible if mtime < deltat
        ! note that currently mtime = deltat
        !================================================================

        if (do_cldice .and. nitend(i,k) > zero .and. ni(i,k)+nitend(i,k)*deltat > nimax(i,k)) then
           nitend(i,k) = max(zero, (nimax(i,k)-ni(i,k))*oneodt)
        end if

     end do

     ! End of "administration" loop

  end do micro_vert_loop ! end k loop

! if (lprnt) write(0,*)' tlat3=',tlat(1,:)*deltat
  !-----------------------------------------------------
  ! convert rain/snow q and N for output to history, note,
  ! output is for gridbox average

  do k=1,nlev
    do i=1,mgncol
      qrout(i,k) = qr(i,k)
      nrout(i,k) = nr(i,k) * rho(i,k)
      qsout(i,k) = qs(i,k)
      nsout(i,k) = ns(i,k) * rho(i,k)
    enddo
  enddo

  ! calculate n0r and lamr from rain mass and number
  ! divide by precip fraction to get in-precip (local) values of
  ! rain mass and number, divide by rhow to get rain number in kg^-1

  do k=1,nlev

     call size_dist_param_basic(mg_rain_props, qric(:,k), nric(:,k), lamr(:,k), mgncol, n0=n0r(:,k))

  enddo
     ! Calculate rercld

     ! calculate mean size of combined rain and cloud water

     call calc_rercld(lamr, n0r, lamc, pgam, qric, qcic, ncic, rercld, mgncol, nlev)


  ! Assign variables back to start-of-timestep values
  ! Some state variables are changed before the main microphysics loop
  ! to make "instantaneous" adjustments. Afterward, we must move those changes
  ! back into the tendencies.
  ! These processes:
  !  - Droplet activation (npccn, impacts nc)
  !  - Instantaneous snow melting  (minstsm/ninstsm, impacts qr/qs/nr/ns)
  !  - Instantaneous rain freezing (minstfr/ninstrf, impacts qr/qs/nr/ns)
  !================================================================================

  do k=1,nlev
    do i=1,mgncol
  ! Re-apply droplet activation tendency
      nc(i,k)     = ncn(i,k)
      nctend(i,k) = nctend(i,k) + npccn(i,k)

  ! Re-apply rain freezing and snow melting.
      qstend(i,k) = qstend(i,k) + (qs(i,k)-qsn(i,k)) * oneodt
      qs(i,k)     = qsn(i,k)

      nstend(i,k) = nstend(i,k) + (ns(i,k)-nsn(i,k)) * oneodt
      ns(i,k)     = nsn(i,k)

      qrtend(i,k) = qrtend(i,k) + (qr(i,k)-qrn(i,k)) * oneodt
      qr(i,k)     = qrn(i,k)

      nrtend(i,k) = nrtend(i,k) + (nr(i,k)-nrn(i,k)) * oneodt
      nr(i,k)     = nrn(i,k)

  !.............................................................................

  !================================================================================

  ! modify to include snow. in prain & evap (diagnostic here: for wet dep)
      nevapr(i,k) = nevapr(i,k) + evapsnow(i,k)
      prain(i,k)  = prain(i,k)  + prodsnow(i,k)

    enddo
  enddo

  do k=1,nlev

     do i=1,mgncol

        ! calculate sedimentation for cloud water and ice
        !================================================================================

        ! update in-cloud cloud mixing ratio and number concentration
        ! with microphysical tendencies to calculate sedimentation, assign to dummy vars
        ! note: these are in-cloud values***, hence we divide by cloud fraction

        if (lcldm(i,k) > mincld) then
          tx1        = one / lcldm(i,k)
          dumc(i,k)  = (qc(i,k)+qctend(i,k)*deltat) * tx1
          dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat)*tx1, zero)
        else
          dumc(i,k)  = zero
          dumnc(i,k) = zero
        endif
        if (icldm(i,k) > mincld) then
          tx1        = one / icldm(i,k)
          dumi(i,k)  = (qi(i,k)+qitend(i,k)*deltat) * tx1
          dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat)*tx1, zero)
        else
          dumi(i,k)  = zero
          dumni(i,k) = zero
        endif
        if (precip_frac(i,k) > mincld) then
          tx1        = one / precip_frac(i,k)
          dumr(i,k)  = (qr(i,k)+qrtend(i,k)*deltat) * tx1
          dums(i,k)  = (qs(i,k)+qstend(i,k)*deltat) * tx1

          dumnr(i,k) = max((nr(i,k)+nrtend(i,k)*deltat)*tx1, zero)
          dumns(i,k) = max((ns(i,k)+nstend(i,k)*deltat)*tx1, zero)
        else
          dumr(i,k)  = zero
          dumr(i,k)  = zero
          dums(i,k)  = zero
          dumns(i,k) = zero
        endif

        ! switch for specification of droplet and crystal number
        if (nccons) then
           dumnc(i,k) = ncnst*rhoinv(i,k)
        end if

        ! switch for specification of cloud ice number
        if (nicons) then
           dumni(i,k) = ninst*rhoinv(i,k)
        end if
     enddo
  enddo

  do k=1,nlev

! obtain new slope parameter to avoid possible singularity


     call size_dist_param_ice(mg_ice_props, dumi(:,k), dumni(:,k), &
                              lami(:,k), mgncol)

     call size_dist_param_liq(mg_liq_props, dumc(:,k), dumnc(:,k), rho(:,k), &
                              pgam(:,k), lamc(:,k), mgncol)

!    call size_dist_param_basic(mg_ice_props, dumi(:,k), dumni(:,k), &
!                               lami(:,k), mgncol)
! fallspeed for rain

     call size_dist_param_basic(mg_rain_props, dumr(:,k), dumnr(:,k), &
                                lamr(:,k), mgncol)
! fallspeed for snow
     call size_dist_param_basic(mg_snow_props, dums(:,k), dumns(:,k), &
                                lams(:,k), mgncol)

  enddo

  do k=1,nlev
     do i=1,mgncol

        ! calculate number and mass weighted fall velocity for droplets and cloud ice
        !-------------------------------------------------------------------

        grho = g*rho(i,k)

        if (dumc(i,k) >= qsmall) then

           tx1 = lamc(i,k)**bc
           vtrmc(i,k) = acn(i,k)*gamma(four+bc+pgam(i,k))          &
                      / (tx1*gamma(pgam(i,k)+four))

           fc(i,k)  = grho*vtrmc(i,k)

           fnc(i,k) = grho* acn(i,k)*gamma(pgam(i,k)+one+bc) &
                              / (tx1*gamma(pgam(i,k)+one))
        else
           fc(i,k)  = zero
           fnc(i,k) = zero
        end if

        ! calculate number and mass weighted fall velocity for cloud ice

        if (dumi(i,k) >= qsmall) then

           tx3 = one / lami(i,k)
           tx1 = ain(i,k) * tx3**bi
           tx2 = 1.2_r8*rhof(i,k)
           vtrmi(i,k) = min(tx1*gamma_bi_plus4*oneo6, tx2)

           fi(i,k)    = grho * vtrmi(i,k)
           fni(i,k)   = grho * min(tx1*gamma_bi_plus1, tx2)

           ! adjust the ice fall velocity for smaller (r < 20 um) ice
           ! particles (blend over 18-20 um)
           irad = (1.5_r8 * 1e6_r8) * tx3
           ifrac = min(one, max(zero, (irad-18._r8)*half))
 
           if (ifrac < one) then
              tx1 = ajn(i,k) / lami(i,k)**bj
              vtrmi(i,k) = ifrac*vtrmi(i,k) + (one-ifrac) * min(tx1*gamma_bj_plus4*oneo6, tx2)

              fi(i,k)    = grho*vtrmi(i,k)
              fni(i,k)   = ifrac * fni(i,k) + (one-ifrac) * grho * min(tx1*gamma_bj_plus1, tx2)
           end if
        else
           fi(i,k) = zero
           fni(i,k)= zero
        end if


        ! fallspeed for rain

!       if (lamr(i,k) >= qsmall) then
        if (dumr(i,k) >= qsmall) then

           ! 'final' values of number and mass weighted mean fallspeed for rain (m/s)

           tx1 = arn(i,k) / lamr(i,k)**br
           tx2 = 9.1_r8*rhof(i,k)
           umr(i,k) = min(tx1*gamma_br_plus4*oneo6, tx2)
           unr(i,k) = min(tx1*gamma_br_plus1,       tx2)

           fr(i,k)  = grho * umr(i,k)
           fnr(i,k) = grho * unr(i,k)

        else
           fr(i,k)  = zero
           fnr(i,k) = zero
        end if

        ! fallspeed for snow


!       if (lams(i,k) >= qsmall) then
        if (dums(i,k) >= qsmall) then

           ! 'final' values of number and mass weighted mean fallspeed for snow (m/s)
           tx1 = asn(i,k) / lams(i,k)**bs
           tx2 = 1.2_r8*rhof(i,k)
           ums(i,k) = min(tx1*gamma_bs_plus4*oneo6, tx2)
           uns(i,k) = min(tx1*gamma_bs_plus1,       tx2)

           fs(i,k)  = grho * ums(i,k)
           fns(i,k) = grho * uns(i,k)

        else
           fs(i,k)  = zero
           fns(i,k) = zero
        end if

        ! redefine dummy variables - sedimentation is calculated over grid-scale
        ! quantities to ensure conservation

        dumc(i,k)  = qc(i,k) + qctend(i,k)*deltat
        dumi(i,k)  = qi(i,k) + qitend(i,k)*deltat
        dumr(i,k)  = qr(i,k) + qrtend(i,k)*deltat
        dums(i,k)  = qs(i,k) + qstend(i,k)*deltat

        dumnc(i,k) = nc(i,k) + nctend(i,k)*deltat
        dumni(i,k) = ni(i,k) + nitend(i,k)*deltat
        dumnr(i,k) = nr(i,k) + nrtend(i,k)*deltat
        dumns(i,k) = ns(i,k) + nstend(i,k)*deltat

        if (dumc(i,k) < qsmall) dumnc(i,k) = zero
        if (dumi(i,k) < qsmall) dumni(i,k) = zero
        if (dumr(i,k) < qsmall) dumnr(i,k) = zero
        if (dums(i,k) < qsmall) dumns(i,k) = zero

     enddo
  end do       !!! vertical loop

  do k=1,nlev
     do i=1,mgncol
       pdel_inv(i,k) = one / pdel(i,k)
     enddo
  enddo
! if (lprnt) write(0,*)' bef sedimentation dumc=',dumc(i,nlev-10:nlev)

  ! initialize nstep for sedimentation sub-steps

  ! calculate number of split time steps to ensure courant stability criteria
  ! for sedimentation calculations
  !-------------------------------------------------------------------
  do i=1,mgncol
     nlb = nlball(i)
     nstep = 1 + nint(max( maxval( fi(i,nlb:nlev)*pdel_inv(i,nlb:nlev)), &
                           maxval(fni(i,nlb:nlev)*pdel_inv(i,nlb:nlev))) * deltat)

     nstep = min(nstep, nstep_def)


     ! loop over sedimentation sub-time step to ensure stability
     !==============================================================
     if (do_cldice) then
       tx2 = one / nstep
       tx1 = tx2 * deltat
       tx3 = tx2 / g

       do n = 1,nstep

        ! top of model

          k = 1

        ! add fallout terms to microphysical tendencies

          tx5           = dumi(i,k)
          tx7           = pdel_inv(i,k) * tx1
          dumi(i,k)     = tx5 / (one + fi(i,k)*tx7)
          tx6           = (dumi(i,k)-tx5) * oneodt
          qitend(i,k)   = qitend(i,k) + tx6
          tx5           = dumni(i,k)
          dumni(i,k)    = tx5 / (one + fni(i,k)*tx7)
          nitend(i,k)   = nitend(i,k) + (dumni(i,k)-tx5) * oneodt

        ! sedimentation tendency for output
          qisedten(i,k) = qisedten(i,k) + tx6

          falouti(k)    = fi(i,k)  * dumi(i,k)
          faloutni(k)   = fni(i,k) * dumni(i,k)

          iflx(i,k+1)   = iflx(i,k+1) + falouti(k) * tx3   ! Ice flux

          do k = 2,nlev

          ! for cloud liquid and ice, if cloud fraction increases with height
          ! then add flux from above to both vapor and cloud water of current level
          ! this means that flux entering clear portion of cell from above evaporates
          ! instantly

          ! note: this is not an issue with precip, since we assume max overlap

            if (icldm(i,k-1) > mincld) then
              dum1 = max(zero, min(one, icldm(i,k)/icldm(i,k-1)))
            else
              dum1 = one
            endif

            tx5           = dumi(i,k)
            tx7           = pdel_inv(i,k) * tx1
            dum2          = tx7 * dum1
            dumi(i,k)     = (tx5 + falouti(k-1)*dum2) / (one + fi(i,k)*tx7)
            tx6           = (dumi(i,k)-tx5) * oneodt
                                                                   ! add fallout terms to eulerian tendencies
            qitend(i,k)   = qitend(i,k) + tx6
            tx5           = dumni(i,k)
            dumni(i,k)    = (tx5 + faloutni(k-1)*dum2)  / (one + fni(i,k)*tx7)
            nitend(i,k)   = nitend(i,k) + (dumni(i,k)-tx5) * oneodt


            qisedten(i,k) = qisedten(i,k) + tx6                    ! sedimentation tendency for output


            falouti(k)  = fi(i,k)  * dumi(i,k)
            faloutni(k) = fni(i,k) * dumni(i,k)

            dum2 = (one-dum1) * falouti(k-1) * pdel_inv(i,k) * tx2
            qvlat(i,k)    = qvlat(i,k)   + dum2                    ! add terms to evap/sub of cloud ice
            qisevap(i,k)  = qisevap(i,k) + dum2                    ! for output

            tlat(i,k)     = tlat(i,k)    - dum2 * xxls

            iflx(i,k+1)   = iflx(i,k+1) + falouti(k) * tx3         ! Ice flux
          end do

        ! units below are m/s
        ! sedimentation flux at surface is added to precip flux at surface
        ! to get total precip (cloud + precip water) rate

          prect(i) = prect(i) + falouti(nlev) * (tx3*0.001_r8)
          preci(i) = preci(i) + falouti(nlev) * (tx3*0.001_r8)

       end do
     end if

! if (lprnt) write(0,*)' tlat4=',tlat(1,:)*deltat
     ! calculate number of split time steps to ensure courant stability criteria
     ! for sedimentation calculations
     !-------------------------------------------------------------------
     nstep = 1 + nint(max( maxval( fc(i,nlb:nlev)*pdel_inv(i,nlb:nlev)), &
                           maxval(fnc(i,nlb:nlev)*pdel_inv(i,nlb:nlev))) * deltat)

     nstep = min(nstep, nstep_def)

     ! loop over sedimentation sub-time step to ensure stability
     !==============================================================
     tx2 = one / nstep
     tx1 = tx2 * deltat
     tx3 = tx2 / g

     do n = 1,nstep

        ! top of model
        k = 1

        tx5           = dumc(i,k)
        tx7           = pdel_inv(i,k) * tx1
        dumc(i,k)     = tx5 / (one + fc(i,k)*tx7)
        tx6           = (dumc(i,k)-tx5) * oneodt
        qctend(i,k)   = qctend(i,k) + tx6
        tx5           = dumnc(i,k)
        dumnc(i,k)    = tx5 / (one + fnc(i,k)*tx7)
        nctend(i,k)   = nctend(i,k) + (dumnc(i,k)-tx5) * oneodt


        ! sedimentation tendency for output
        qcsedten(i,k) = qcsedten(i,k) + tx6

        faloutc(k)    = fc(i,k)  * dumc(i,k)
        faloutnc(k)   = fnc(i,k) * dumnc(i,k)

        lflx(i,k+1)   = lflx(i,k+1) + faloutc(k) * tx3
        do k = 2,nlev

          if (lcldm(i,k-1) > mincld) then
            dum1 = max(zero, min(one, lcldm(i,k)/lcldm(i,k-1)))
          else
            dum1 = one
          endif

          tx5            = dumc(i,k)
          tx7            = pdel_inv(i,k) * tx1
          dum2           = tx7 * dum1
          dumc(i,k)      = (tx5 + faloutc(k-1)*dum2) / (one + fc(i,k)*tx7)
          tx6            = (dumc(i,k)-tx5) * oneodt
          qctend(i,k)    = qctend(i,k) + tx6
          tx5            = dumnc(i,k)
          dumnc(i,k)     = (tx5 + faloutnc(k-1)*dum2)  / (one + fnc(i,k)*tx7)
          nctend(i,k)    = nctend(i,k) + (dumnc(i,k)-tx5) * oneodt



           qcsedten(i,k) = qcsedten(i,k) + tx6               ! sedimentation tendency for output

           faloutc(k)  = fc(i,k)  * dumc(i,k)
           faloutnc(k) = fnc(i,k) * dumnc(i,k)

           dum2 = (one-dum1) * faloutc(k-1) * pdel_inv(i,k) * tx2
           qvlat(i,k)    = qvlat(i,k)   + dum2               ! add terms to to evap/sub of cloud water
           qcsevap(i,k)  = qcsevap(i,k) + dum2               ! for output

           tlat(i,k)     = tlat(i,k)    - dum2 * xxlv

           lflx(i,k+1)   = lflx(i,k+1) + faloutc(k) * tx3    ! Liquid condensate flux here
        end do

        prect(i) = prect(i) + faloutc(nlev) * (tx3*0.001_r8)

     end do
! if (lprnt) write(0,*)' tlat5=',tlat(1,:)*deltat
! if (lprnt) write(0,*)' maxval=',maxval( fr(i,nlb:nlev)*pdel_inv(i,nlb:nlev))&
!                                ,maxval(fnr(i,nlb:nlev)*pdel_inv(i,nlb:nlev))

     ! calculate number of split time steps to ensure courant stability criteria
     ! for sedimentation calculations
     !-------------------------------------------------------------------
     nstep = 1 + nint(max( maxval( fr(i,nlb:nlev)*pdel_inv(i,nlb:nlev)), &
                           maxval(fnr(i,nlb:nlev)*pdel_inv(i,nlb:nlev))) * deltat)

     nstep = min(nstep, nstep_def)

     ! loop over sedimentation sub-time step to ensure stability
     !==============================================================
     tx2 = one / nstep
     tx1 = tx2 * deltat
     tx3 = tx2 / g

!    if(lprnt) then
!      write(0,*)' nstep=',nstep,' tx1=',tx1,' tx2=',tx2,' tx3=',tx3,' qsmall=',qsmall
!      write(0,*)' fr=',fr(i,:)
!      write(0,*)' dumr=',dumr(i,:)
!    endif

     do n = 1,nstep

        ! top of model
        k = 1

        ! add fallout terms to microphysical tendencies

        tx5 = dumr(i,k)
        tx7 = pdel_inv(i,k) * tx1
        dumr(i,k)   = tx5 / (one + fr(i,k)*tx7)
        tx6 = (dumr(i,k)-tx5) * oneodt
        qrtend(i,k) = qrtend(i,k) + tx6
        tx5 = dumnr(i,k)
        dumnr(i,k)  = tx5 / (one + fnr(i,k)*tx7)
        nrtend(i,k) = nrtend(i,k) + (dumnr(i,k)-tx5) * oneodt

        ! sedimentation tendency for output
        qrsedten(i,k) = qrsedten(i,k) + tx6

        faloutr(k)  = fr(i,k)  * dumr(i,k)
        faloutnr(k) = fnr(i,k) * dumnr(i,k)

        rflx(i,k+1) = rflx(i,k+1) + faloutr(k) * tx3

        do k = 2,nlev

          tx5 = dumr(i,k)
          tx7 = pdel_inv(i,k) * tx1
          dumr(i,k)   = (tx5 + faloutr(k-1)*tx7) / (one + fr(i,k)*tx7)
          tx6 = (dumr(i,k)-tx5) * oneodt
          qrtend(i,k) = qrtend(i,k) + tx6
          tx5 = dumnr(i,k)
          dumnr(i,k)  = (tx5 + faloutnr(k-1)*tx7)  / (one + fnr(i,k)*tx7)
          nrtend(i,k) = nrtend(i,k) + (dumnr(i,k)-tx5) * oneodt


                                                          ! sedimentation tendency for output
           qrsedten(i,k) = qrsedten(i,k) + tx6                      ! sedimentation tendency for output

           faloutr(k)  = fr(i,k)  * dumr(i,k)
           faloutnr(k) = fnr(i,k) * dumnr(i,k)

           rflx(i,k+1)   = rflx(i,k+1) + faloutr(k) * tx3 ! Rain Flux
        end do

        prect(i) = prect(i) + faloutr(nlev) * (tx3*0.001_r8)

     end do

! if (lprnt) write(0,*)' prectaftrain=',prect(i),' preci=',preci(i)
     ! calculate number of split time steps to ensure courant stability criteria
     ! for sedimentation calculations
     !-------------------------------------------------------------------
     nstep = 1 + nint(max( maxval( fs(i,nlb:nlev)*pdel_inv(i,nlb:nlev)), &
                           maxval(fns(i,nlb:nlev)*pdel_inv(i,nlb:nlev))) * deltat)
     nstep = min(nstep, nstep_def)

     ! loop over sedimentation sub-time step to ensure stability
     !==============================================================
     tx2 = one / nstep
     tx1 = tx2 * deltat
     tx3 = tx2 / g
     do n = 1,nstep

        ! top of model
        k = 1

        ! add fallout terms to microphysical tendencies

        tx5 = dums(i,k)
        tx7 = pdel_inv(i,k) * tx1
        dums(i,k)   = tx5 / (one + fs(i,k)*tx7)
        tx6 = (dums(i,k)-tx5) * oneodt
        qstend(i,k) = qstend(i,k) + tx6
        tx5 = dumns(i,k)
        dumns(i,k)  = tx5 / (one + fns(i,k)*tx7)
        nstend(i,k) = nstend(i,k) + (dumns(i,k)-tx5) * oneodt

        ! sedimentation tendency for output
        qssedten(i,k) = qssedten(i,k) + tx6

        falouts(k)  = fs(i,k)  * dums(i,k)
        faloutns(k) = fns(i,k) * dumns(i,k)

        sflx(i,k+1)   = sflx(i,k+1) + falouts(k) * tx3

        do k = 2,nlev


           tx5 = dums(i,k)
           tx7 = pdel_inv(i,k) * tx1
           dums(i,k)   = (tx5 + falouts(k-1)*tx7) / (one + fs(i,k)*tx7)
           tx6 = (dums(i,k)-tx5) * oneodt
           qstend(i,k) = qstend(i,k) + tx6
           tx5 = dumns(i,k)
           dumns(i,k)  = (tx5 + faloutns(k-1)*tx7)  / (one + fns(i,k)*tx7)
           nstend(i,k) = nstend(i,k) + (dumns(i,k)-tx5) * oneodt


           qssedten(i,k) = qssedten(i,k) + tx6                      ! sedimentation tendency for output

           falouts(k)  = fs(i,k)  * dums(i,k)
           faloutns(k) = fns(i,k) * dumns(i,k)

           sflx(i,k+1)   = sflx(i,k+1) + falouts(k) * tx3 ! Snow Flux
        end do   !! k loop

        prect(i) = prect(i) + falouts(nlev) * (tx3*0.001_r8)
        preci(i) = preci(i) + falouts(nlev) * (tx3*0.001_r8)

     end do   !! nstep loop

  enddo       ! end of i loop
  ! end sedimentation

! if (lprnt) write(0,*)' prectaftsed=',prect(i),' preci=',preci(i)

  !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

  ! get new update for variables that includes sedimentation tendency
  ! note : here dum variables are grid-average, NOT in-cloud

  do k=1,nlev
     do i=1,mgncol
        dumc(i,k)  = max(qc(i,k)+qctend(i,k)*deltat, zero)
        dumi(i,k)  = max(qi(i,k)+qitend(i,k)*deltat, zero)
        dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat, zero)
        dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat, zero)

        dumr(i,k)  = max(qr(i,k)+qrtend(i,k)*deltat, zero)
        dumnr(i,k) = max(nr(i,k)+nrtend(i,k)*deltat, zero)
        dums(i,k)  = max(qs(i,k)+qstend(i,k)*deltat, zero)
        dumns(i,k) = max(ns(i,k)+nstend(i,k)*deltat, zero)

        ! switch for specification of droplet and crystal number
        if (nccons) then
           dumnc(i,k) = ncnst*rhoinv(i,k)*lcldm(i,k)
        end if

        ! switch for specification of cloud ice number
        if (nicons) then
           dumni(i,k) = ninst*rhoinv(i,k)*icldm(i,k)
        end if

        if (dumc(i,k) < qsmall) dumnc(i,k) = zero
        if (dumi(i,k) < qsmall) dumni(i,k) = zero
        if (dumr(i,k) < qsmall) dumnr(i,k) = zero
        if (dums(i,k) < qsmall) dumns(i,k) = zero

     enddo

  enddo

  ! calculate instantaneous processes (melting, homogeneous freezing)
  !====================================================================

  ! melting of snow at +2 C
  do k=1,nlev

     do i=1,mgncol

        tx1 = t(i,k) + tlat(i,k)*(deltat/cpp) - snowmelt
        if (tx1 > zero) then
           if (dums(i,k) > zero) then

              ! make sure melting snow doesn't reduce temperature below threshold
              dum = -(xlf/cpp) * dums(i,k)
              if (tx1+dum <  zero) then
                 dum = min(one, max(zero, -tx1/dum))
              else
                 dum = one
              end if

              tx1 = dum * oneodt
              qstend(i,k) = qstend(i,k) - tx1*dums(i,k)
              nstend(i,k) = nstend(i,k) - tx1*dumns(i,k)
              qrtend(i,k) = qrtend(i,k) + tx1*dums(i,k)
              nrtend(i,k) = nrtend(i,k) + tx1*dumns(i,k)

              dum1 = - xlf * tx1 * dums(i,k)
              tlat(i,k)       = tlat(i,k)       + dum1
              meltsdttot(i,k) = meltsdttot(i,k) + dum1
           end if
        end if
     enddo
  enddo
   do k=1,nlev
      do i=1,mgncol

        ! freezing of rain at -5 C

        tx1 = t(i,k) + tlat(i,k) * (deltat/cpp) - rainfrze
        if (tx1 < zero) then

           if (dumr(i,k) > zero) then

              ! make sure freezing rain doesn't increase temperature above threshold
              dum = (xlf/cpp) * dumr(i,k)
              if (tx1+dum > zero) then
                 dum = min(one, max(zero, -tx1/dum))
              else
                 dum = one
              end if
              tx2 = dum * oneodt
              qrtend(i,k) = qrtend(i,k) - tx2 * dumr(i,k)
              nrtend(i,k) = nrtend(i,k) - tx2 * dumnr(i,k)

              ! get mean size of rain = 1/lamr, add frozen rain to either snow or cloud ice
              ! depending on mean rain size

              call size_dist_param_basic(mg_rain_props, dumr(i,k), dumnr(i,k), lamr(i,k))

              if (lamr(i,k) < one/Dcs) then
                 qstend(i,k) = qstend(i,k) + tx2 * dumr(i,k)
                 nstend(i,k) = nstend(i,k) + tx2 * dumnr(i,k)
              else
                 qitend(i,k) = qitend(i,k) + tx2 * dumr(i,k)
                 nitend(i,k) = nitend(i,k) + tx2 * dumnr(i,k)
              end if
                                                         ! heating tendency
              dum1 = xlf*dum*dumr(i,k)*oneodt
              frzrdttot(i,k) = dum1 + frzrdttot(i,k)
              tlat(i,k)      = dum1 + tlat(i,k)

           end if
        end if

      enddo
   enddo
   if (do_cldice) then
      do k=1,nlev
        do i=1,mgncol
           tx1 = t(i,k) + tlat(i,k) * (deltat/cpp) - tmelt
           if (tx1 > zero) then
              if (dumi(i,k) > zero) then

                 ! limit so that melting does not push temperature below freezing
                 !-----------------------------------------------------------------
                 dum = -dumi(i,k)*xlf/cpp
                 if (tx1+dum < zero) then
                    dum = min(one, max(zero, -tx1/dum))
                 else
                    dum = one
                 end if

                 tx2 = dum * oneodt
                 qctend(i,k) = qctend(i,k) + tx2*dumi(i,k)

                 ! for output
                 melttot(i,k) = tx2*dumi(i,k)

                 ! assume melting ice produces droplet
                 ! mean volume radius of 8 micron

                 nctend(i,k) = nctend(i,k) + three*tx2*dumi(i,k)/(four*pi*5.12e-16_r8*rhow)

                 qitend(i,k) = ((one-dum)*dumi(i,k)-qi(i,k))  * oneodt
                 nitend(i,k) = ((one-dum)*dumni(i,k)-ni(i,k)) * oneodt
                 tlat(i,k)   = tlat(i,k) - xlf*tx2*dumi(i,k)
              end if
           end if
        enddo
     enddo

! if (lprnt) write(0,*)' tlat6=',tlat(1,:)*deltat
! if (lprnt) write(0,*)' qitend=',qitend(1,nlev-10:nlev)*deltat
! if (lprnt) write(0,*)' qctend=',qctend(1,nlev-10:nlev)*deltat
     ! homogeneously freeze droplets at -40 C
     !-----------------------------------------------------------------

     do k=1,nlev
        do i=1,mgncol
           tx1 = t(i,k) + tlat(i,k)*(deltat/cpp) - 233.15_r8
           if (tx1 < zero) then
              if (dumc(i,k) > zero) then

                 ! limit so that freezing does not push temperature above threshold
                 dum = (xlf/cpp) * dumc(i,k)
                 if (tx1+dum > zero) then
                    dum = min(one, max(zero, -tx1/dum))
                 else
                    dum = one
                 end if

                 tx2 = dum * oneodt * dumc(i,k)
                 qitend(i,k)  = tx2 + qitend(i,k)
                 homotot(i,k) = tx2                   ! for output

                 ! assume 25 micron mean volume radius of homogeneously frozen droplets
                 ! consistent with size of detrained ice in stratiform.F90

                 nitend(i,k) = nitend(i,k) + tx2*(three/(four*pi*1.563e-14_r8* 500._r8))
                 qctend(i,k) = ((one-dum)*dumc(i,k)-qc(i,k))  * oneodt
                 nctend(i,k) = ((one-dum)*dumnc(i,k)-nc(i,k)) * oneodt
                 tlat(i,k)   = tlat(i,k) + xlf*tx2
              end if
           end if
        enddo 
     enddo 
     ! remove any excess over-saturation, which is possible due to non-linearity when adding
     ! together all microphysical processes
     !-----------------------------------------------------------------
     ! follow code similar to old CAM scheme
     do k=1,nlev
        do i=1,mgncol

           qtmp = q(i,k) + qvlat(i,k) * deltat
           ttmp = t(i,k) + tlat(i,k)  * (deltat/cpp)

           ! use rhw to allow ice supersaturation
           !call qsat_water(ttmp, p(i,k), esn, qvn)
           esn = min(fpvsl(ttmp), p(i,k))
           qvn = epsqs*esn/(p(i,k)-omeps*esn) * qsfm(i,k)


           if (qtmp > qvn .and. qvn > 0 .and. allow_sed_supersat) then
              ! expression below is approximate since there may be ice deposition
              dum = (qtmp-qvn)/(one+xxlv_squared*qvn/(cpp*rv*ttmp*ttmp)) * oneodt
              ! add to output cme
              cmeout(i,k) = cmeout(i,k) + dum
              ! now add to tendencies, partition between liquid and ice based on temperature
              if (ttmp > 268.15_r8) then
                 dum1 = zero
                 ! now add to tendencies, partition between liquid and ice based on te
                 !-------------------------------------------------------
              else if (ttmp < 238.15_r8) then
                 dum1 = one
              else
                 dum1 = (268.15_r8-ttmp)/30._r8
              end if

              tx1 = xxls*dum1 + xxlv*(one-dum1)
              dum = (qtmp-qvn)/(one+tx1*tx1*qvn/(cpp*rv*ttmp*ttmp)) * oneodt
              tx2 = dum*(one-dum1)
              qctend(i,k)   = qctend(i,k) + tx2
              qcrestot(i,k) = tx2                         ! for output
              qitend(i,k)   = qitend(i,k) + dum*dum1
              qirestot(i,k) = dum*dum1
              qvlat(i,k)    = qvlat(i,k) - dum
              ! for output
              qvres(i,k)    = -dum
              tlat(i,k)     = tlat(i,k) + dum*tx1
           end if
        enddo 
     enddo 
  end if

! if (lprnt) write(0,*)' tlat7=',tlat(1,:)*deltat
  ! calculate effective radius for pass to radiation code
  !=========================================================
  ! if no cloud water, default value is 10 micron for droplets,
  ! 25 micron for cloud ice

  ! update cloud variables after instantaneous processes to get effective radius
  ! variables are in-cloud to calculate size dist parameters
  do k=1,nlev
     do i=1,mgncol
        if (lcldm(i,k) > mincld) then
          tx1 = one / lcldm(i,k)
        else
          tx1 = zero
        endif
        if (icldm(i,k) > mincld) then
          tx2 = one / icldm(i,k)
        else
          tx2 = zero
        endif
        if (precip_frac(i,k) > mincld) then
          tx3 = one / precip_frac(i,k)
        else
          tx3 = zero
        endif
        dumc(i,k)  = max(qc(i,k)+qctend(i,k)*deltat, zero) * tx1
        dumi(i,k)  = max(qi(i,k)+qitend(i,k)*deltat, zero) * tx2
        dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat, zero) * tx1
        dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat, zero) * tx2

        dumr(i,k)  = max(qr(i,k)+qrtend(i,k)*deltat, zero) * tx3
        dumnr(i,k) = max(nr(i,k)+nrtend(i,k)*deltat, zero) * tx3
        dums(i,k)  = max(qs(i,k)+qstend(i,k)*deltat, zero) * tx3
        dumns(i,k) = max(ns(i,k)+nstend(i,k)*deltat, zero) * tx3

        ! switch for specification of droplet and crystal number
        if (nccons) then
           dumnc(i,k) = ncnst * rhoinv(i,k)
        end if

        ! switch for specification of cloud ice number
        if (nicons) then
           dumni(i,k) = ninst * rhoinv(i,k)
        end if

        ! limit in-cloud mixing ratio to reasonable value of 5 g kg-1
!       dumc(i,k) = min(dumc(i,k), 5.e-3_r8)
!       dumi(i,k) = min(dumi(i,k), 5.e-3_r8)
        dumc(i,k) = min(dumc(i,k), 10.e-3_r8)
        dumi(i,k) = min(dumi(i,k), 10.e-3_r8)
        ! limit in-precip mixing ratios
        dumr(i,k) = min(dumr(i,k), 10.e-3_r8)
        dums(i,k) = min(dums(i,k), 10.e-3_r8)
     enddo
  enddo
  ! cloud ice effective radius
  !-----------------------------------------------------------------

  if (do_cldice) then
     do k=1,nlev
        do i=1,mgncol
           if (dumi(i,k) >= qsmall) then

              tx1 = dumni(i,k)
              call size_dist_param_basic(mg_ice_props, dumi(i,k), dumni(i,k), &
                                         lami(i,k), dumni0)

              if (dumni(i,k) /= tx1) then
                 ! adjust number conc if needed to keep mean size in reasonable range
                 nitend(i,k) = (dumni(i,k)*icldm(i,k)-ni(i,k)) * oneodt
              end if

              tx1 = one / lami(i,k)
!             effi(i,k)   = (1.5_r8*1.e6_r8) * tx1
              effi(i,k)   = (three*1.e6_r8) * tx1
              sadice(i,k) = two*pi*(tx1*tx1*tx1)*dumni0*rho(i,k)*1.e-2_r8  ! m2/m3 -> cm2/cm3

           else
!             effi(i,k)   = 25._r8
              effi(i,k)   = 50._r8
              sadice(i,k) = zero
           end if

           ! ice effective diameter for david mitchell's optics
           deffi(i,k) = effi(i,k) * (rhoi+rhoi)/rhows
        enddo
     enddo
  !else
     !do k=1,nlev
        !do i=1,mgncol
           ! NOTE: If CARMA is doing the ice microphysics, then the ice effective
           ! radius has already been determined from the size distribution.
           !effi(i,k) = re_ice(i,k) * 1.e6_r8      ! m -> um
           !deffi(i,k)=effi(i,k) * 2._r8
           !sadice(i,k) = 4._r8*pi*(effi(i,k)**2)*ni(i,k)*rho(i,k)*1e-2_r8
        !enddo
     !enddo
  end if

  ! cloud droplet effective radius
  !-----------------------------------------------------------------
  do k=1,nlev
     do i=1,mgncol
        if (dumc(i,k) >= qsmall) then


           ! switch for specification of droplet and crystal number
           if (nccons) then
              ! make sure nc is consistence with the constant N by adjusting tendency, need
              ! to multiply by cloud fraction
              ! note that nctend may be further adjusted below if mean droplet size is
              ! out of bounds

              nctend(i,k) = (ncnst*rhoinv(i,k)*lcldm(i,k)-nc(i,k)) * oneodt

           end if

           dum = dumnc(i,k)

           call size_dist_param_liq(mg_liq_props, dumc(i,k), dumnc(i,k), rho(i,k), &
                                    pgam(i,k), lamc(i,k))

           if (dum /= dumnc(i,k)) then
              ! adjust number conc if needed to keep mean size in reasonable range
              nctend(i,k) = (dumnc(i,k)*lcldm(i,k)-nc(i,k)) * oneodt
           end if

           effc(i,k) = (half*1.e6_r8) * (pgam(i,k)+three) / lamc(i,k)
           !assign output fields for shape here
           lamcrad(i,k) = lamc(i,k)
           pgamrad(i,k) = pgam(i,k)


           ! recalculate effective radius for constant number, in order to separate
           ! first and second indirect effects
           !======================================
           ! assume constant number of 10^8 kg-1

           dumnc(i,k) = 1.e8_r8

           ! Pass in "false" adjust flag to prevent number from being changed within
           ! size distribution subroutine.
           call size_dist_param_liq(mg_liq_props, dumc(i,k), dumnc(i,k), rho(i,k), &
                                    pgam(i,k), lamc(i,k))

           effc_fn(i,k) = (half*1.e6_r8) * (pgam(i,k)+three)/lamc(i,k)

        else
           effc(i,k)    = ten
           lamcrad(i,k) = zero
           pgamrad(i,k) = zero
           effc_fn(i,k) = ten
        end if
     enddo
  enddo
  ! recalculate 'final' rain size distribution parameters
  ! to ensure that rain size is in bounds, adjust rain number if needed
  do k=1,nlev
     do i=1,mgncol

        if (dumr(i,k) >= qsmall) then

           dum = dumnr(i,k)

           call size_dist_param_basic(mg_rain_props, dumr(i,k), dumnr(i,k), lamr(i,k))

           if (dum /= dumnr(i,k)) then
              ! adjust number conc if needed to keep mean size in reasonable range
              nrtend(i,k) = (dumnr(i,k)*precip_frac(i,k)-nr(i,k)) *oneodt
           end if

        end if
     enddo
  enddo
  ! recalculate 'final' snow size distribution parameters
  ! to ensure that snow size is in bounds, adjust snow number if needed
  do k=1,nlev
     do i=1,mgncol
        if (dums(i,k) >= qsmall) then

           dum = dumns(i,k)

           call size_dist_param_basic(mg_snow_props, dums(i,k), dumns(i,k), &
                                      lams(i,k), n0=dumns0)

           if (dum /= dumns(i,k)) then
              ! adjust number conc if needed to keep mean size in reasonable range
              nstend(i,k) = (dumns(i,k)*precip_frac(i,k)-ns(i,k)) * oneodt
           end if

           tx1 = (two*pi*1.e-2_r8) / (lams(i,k)*lams(i,k)*lams(i,k))
           sadsnow(i,k) = tx1*dumns0*rho(i,k)             ! m2/m3 -> cm2/cm3

        end if


     end do ! vertical k loop
  enddo
  do k=1,nlev
     do i=1,mgncol
        ! if updated q (after microphysics) is zero, then ensure updated n is also zero
        !=================================================================================
        if (qc(i,k)+qctend(i,k)*deltat < qsmall)                 nctend(i,k) = -nc(i,k) * oneodt
        if (do_cldice .and. qi(i,k)+qitend(i,k)*deltat < qsmall) nitend(i,k) = -ni(i,k) * oneodt
        if (qr(i,k)+qrtend(i,k)*deltat < qsmall)                 nrtend(i,k) = -nr(i,k) * oneodt
        if (qs(i,k)+qstend(i,k)*deltat < qsmall)                 nstend(i,k) = -ns(i,k) * oneodt

     end do

  end do

  ! DO STUFF FOR OUTPUT:
  !==================================================

  do k=1,nlev
    do i=1,mgncol

  ! qc and qi are only used for output calculations past here,
  ! so add qctend and qitend back in one more time
      qc(i,k) = qc(i,k) + qctend(i,k)*deltat
      qi(i,k) = qi(i,k) + qitend(i,k)*deltat

  ! averaging for snow and rain number and diameter
  !--------------------------------------------------

  ! drout2/dsout2:
  ! diameter of rain and snow
  ! dsout:
  ! scaled diameter of snow (passed to radiation in CAM)
  ! reff_rain/reff_snow:
  ! calculate effective radius of rain and snow in microns for COSP using Eq. 9 of COSP v1.3 manual

      if (qrout(i,k) > 1.e-7_r8 .and. nrout(i,k) > zero) then
         qrout2(i,k)    = qrout(i,k) * precip_frac(i,k)
         nrout2(i,k)    = nrout(i,k) * precip_frac(i,k)
     ! The avg_diameter call does the actual calculation; other diameter
     ! outputs are just drout2 times constants.
         drout2(i,k)    = avg_diameter(qrout(i,k), nrout(i,k), rho(i,k), rhow)
         freqr(i,k)     = precip_frac(i,k)

         reff_rain(i,k) = (1.e6_r8*three) * drout2(i,k)
      else
         qrout2(i,k)    = zero
         nrout2(i,k)    = zero
         drout2(i,k)    = zero
         freqr(i,k)     = zero
         reff_rain(i,k) = zero
      endif

      if (qsout(i,k) > 1.e-7_r8 .and. nsout(i,k) > zero) then
         qsout2(i,k)    = qsout(i,k) * precip_frac(i,k)
         nsout2(i,k)    = nsout(i,k) * precip_frac(i,k)
     ! The avg_diameter call does the actual calculation; other diameter
     ! outputs are just dsout2 times constants.
         dsout2(i,k)    = avg_diameter(qsout(i,k), nsout(i,k), rho(i,k), rhosn)
         freqs(i,k)     = precip_frac(i,k)

         dsout(i,k)     = three*rhosn/rhows*dsout2(i,k)

         reff_snow(i,k) = (1.e6_r8*three) * dsout2(i,k)
      else
         dsout(i,k)     = zero
         qsout2(i,k)    = zero
         nsout2(i,k)    = zero 
         dsout2(i,k)    = zero
         freqs(i,k)     = zero
         reff_snow(i,k) = zero
      endif

    enddo
  enddo

  ! analytic radar reflectivity
  !--------------------------------------------------
  ! formulas from Matthew Shupe, NOAA/CERES
  ! *****note: radar reflectivity is local (in-precip average)
  ! units of mm^6/m^3

  do k=1,nlev
     do i = 1,mgncol
        if (qc(i,k) >= qsmall .and. (nc(i,k)+nctend(i,k)*deltat) > ten) then
           tx1 = rho(i,k) / lcldm(i,k)
           tx2 = 1000._r8 * qc(i,k) * tx1
           dum = tx2 * tx2 * lcldm(i,k)                            &
                /(0.109_r8*(nc(i,k)+nctend(i,k)*deltat)*tx1*1.e-6_r8*precip_frac(i,k))
!          dum = (qc(i,k)/lcldm(i,k)*rho(i,k)*1000._r8)**2 &
!               /(0.109_r8*(nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k)*rho(i,k)/1.e6_r8)*lcldm(i,k)/precip_frac(i,k)
        else
           dum = zero
        end if
        if (qi(i,k) >= qsmall) then
!          dum1 = (qi(i,k)*rho(i,k)/icldm(i,k)*1000._r8/0.1_r8)**(one/0.63_r8)*icldm(i,k)/precip_frac(i,k)
           dum1 = (qi(i,k)*rho(i,k)/icldm(i,k)*10000._r8)**(one/0.63_r8)*icldm(i,k)/precip_frac(i,k)
        else
           dum1 = zero
        end if

        if (qsout(i,k) >= qsmall) then
!          dum1 = dum1 + (qsout(i,k)*rho(i,k)*1000._r8/0.1_r8)**(one/0.63_r8)
           dum1 = dum1 + (qsout(i,k)*rho(i,k)*10000._r8)**(one/0.63_r8)
        end if

        refl(i,k) = dum + dum1

        ! add rain rate, but for 37 GHz formulation instead of 94 GHz
        ! formula approximated from data of Matrasov (2007)
        ! rainrt is the rain rate in mm/hr
        ! reflectivity (dum) is in DBz

        if (rainrt(i,k) >= 0.001_r8) then
           dum = rainrt(i,k) * rainrt(i,k)
           dum = log10(dum*dum*dum) + 16._r8

           ! convert from DBz to mm^6/m^3

           dum = ten**(dum/ten)
        else
           ! don't include rain rate in R calculation for values less than 0.001 mm/hr
           dum = zero
        end if

        ! add to refl

        refl(i,k) = refl(i,k) + dum

        !output reflectivity in Z.
        areflz(i,k) = refl(i,k) * precip_frac(i,k)

        ! convert back to DBz

        if (refl(i,k) > minrefl) then
           refl(i,k) = ten*log10(refl(i,k))
        else
           refl(i,k) = -9999._r8
        end if

        !set averaging flag
        if (refl(i,k) > mindbz) then
           arefl(i,k) = refl(i,k) * precip_frac(i,k)
           frefl(i,k) = precip_frac(i,k)
        else
           arefl(i,k)  = zero
           areflz(i,k) = zero
           frefl(i,k)  = zero
        end if

        ! bound cloudsat reflectivity

        csrfl(i,k) = min(csmax,refl(i,k))

        !set averaging flag
        if (csrfl(i,k) > csmin) then
           acsrfl(i,k) = refl(i,k) * precip_frac(i,k)
           fcsrfl(i,k) = precip_frac(i,k)
        else
           acsrfl(i,k) = zero
           fcsrfl(i,k) = zero
        end if

     end do
  end do

  do k=1,nlev
     do i = 1,mgncol
  !redefine fice here....
       tx2 = qsout(i,k) + qi(i,k)
       tx1 = tx2 + qrout(i,k) + qc(i,k)
       if ( tx2 > qsmall .and. tx1 > qsmall) then
         nfice(i,k) = min(tx2/tx1, one)
       else
         nfice(i,k) = zero
       endif
    enddo
  enddo

end subroutine micro_mg_tend

!========================================================================
!OUTPUT CALCULATIONS
!========================================================================

!>\ingroup mg2_0_mp
!! This subroutine 
subroutine calc_rercld(lamr, n0r, lamc, pgam, qric, qcic, ncic, rercld, mgncol,nlev)
  integer, intent(in) :: mgncol, nlev
  real(r8), dimension(mgncol,nlev), intent(in) :: lamr          ! rain size parameter (slope)
  real(r8), dimension(mgncol,nlev), intent(in) :: n0r           ! rain size parameter (intercept)
  real(r8), dimension(mgncol,nlev), intent(in) :: lamc          ! size distribution parameter (slope)
  real(r8), dimension(mgncol,nlev), intent(in) :: pgam          ! droplet size parameter
  real(r8), dimension(mgncol,nlev), intent(in) :: qric          ! in-cloud rain mass mixing ratio
  real(r8), dimension(mgncol,nlev), intent(in) :: qcic          ! in-cloud cloud liquid
  real(r8), dimension(mgncol,nlev), intent(in) :: ncic          ! in-cloud droplet number concentration

  real(r8), dimension(mgncol,nlev), intent(inout) :: rercld     ! effective radius calculation for rain + cloud

  ! combined size of precip & cloud drops
  real(r8) :: Atmp

  integer :: i, k

  do k=1,nlev
    do i=1,mgncol
      ! Rain drops
      if (lamr(i,k) > zero) then
        Atmp = n0r(i,k) * (half*pi) / (lamr(i,k)*lamr(i,k)*lamr(i,k))
      else
        Atmp = zero
      end if

      ! Add cloud drops
      if (lamc(i,k) > zero) then
        Atmp = Atmp + ncic(i,k) * pi * rising_factorial(pgam(i,k)+one, 2) &
                    / (four*lamc(i,k)*lamc(i,k))
      end if

      if (Atmp > zero) then
        rercld(i,k) = rercld(i,k) + three *(qric(i,k) + qcic(i,k)) / (four * rhow * Atmp)
      end if
    enddo
  enddo
end subroutine calc_rercld

!========================================================================

end module micro_mg2_0
